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1.  INTRODUCTION 


For  the  past  seven  years,  A.R.A.P.,  Inc.  has  had  a  series  of  contracts 
from  the  Naval  Air  Systems  Command  to  develop  a  computer  model  for  determining 
the  detailed  low-level  atmospheric  distributions  of  velocity,  temperature, 
moisture,  refractive  index,  and  the  turbulent  variances  of  these  quantities 
for  marine  environments.  In  addition  to  appropriately  modeling  the  turbulent 
transport  of  momentum,  heat,  and  moisture,  it  was  necessary  to  incorporate 
moisture  change  of  phase  and  the  physics  of  thermal  radiation  into  this  model 
since  low-level  clouds  or  fog  are  a  frequent  occurance  in  the  marine 
atmospheric  boundary  layer.  The  development  and  a  number  of  sample 
calculations  exemplifying  different  phenomena  are  detailed  in  References  1-17. 
Reference  18  which  accompanies  this  report  provides  a  detailed  review  of 
modeling  the  atmospheric  boundary  layer  using  turbulent  transport  theory.  It 
includes  a  review  of  the  status  of  our  understanding  of  atmospheric  boundary 
layer  dynamics,  as  well  as  a  review  of  the  modeling  of  the  three  physical 
processes  most  critical  for  determining  the  atmospheric  marine  boundary  layer. 
These  three  are  turbulent  transport,  thermal  radiation,  and  change  of  phase  of 
atmospheric  water.  Reference  18  also  provides  a  review  of  many  of  the  sample 
calculations  made  with  the  A.R.A.P.  model  which  successfully  illustrates 
features  expected  in  the  atmospheric  marine  boundary  layer.  Section  5  of  the 
present  report  counterbalances  those  "successes"  by  detailing  some  problem 
areas  of  the  current  model.  Together  they  provide  the  detailed  critical 
review  specified  in  the  past  years  contract. 

Two  model  calculations  performed  during  the  past  year  are  detailed  in 
Section  2  and  in  Appendix  A.  The  fog  calculation  presented  in  the  next 
section  gives  one  possible  reason  why  the  surface  air  in  fogs  is  generally 
found  to  be  cooler  than  the  ocean  surface.  By  following  the  evolution  of  a 
fog  which  is  formed  by  warm  air  passing  over  colder  water  under  nocturnal 
radiation  conditions,  it  is  shown  that  the  enhanced  radiational  cooling 
induced  by  the  fog  is  sufficient  to  reduce  the  surface  air  temperature 
relatively  rapidly  to  below  that  of  the  water  surface.  During  the  past  year 
we  have  had  discussions  with  E.  Mack  and  W.  Rodgers  of  Cal  span  regarding  the 
forthcoming  fog  model  evaluation  study.  We  look  forward  to  exercising  our 
model  as  part  of  this  study. 

Appendix  A  details  a  calculation  of  the  detailed  mechanism  of 
Kel vin-Helmholtz  wave  breaking.  The  turbulent  breaking  process  is  modeled 
using  our  second-order  closure  model  to  describe  the  small-scale  turbulence, 
while  the  large  scale  billow  itself  is  calculated  explicitly  as  a 
two-dimensional  flow.  This  calculation  was  partially  supported  by  an  ONR 
contract  which  cabled  for  examining  trackable  clear-air,  radar  signals.  The 
large  values  of  C^,  which  can  occur  when  this  phenomenon  occurs  at  the  top  of 
a  relatively  moist  boundary  layer  when  the  air  above  the  inversion  is  much 
dryer,  makes  this  a  likely  candidate.  This  type  of  calculation  can  also  be 
used  to  investigate  the  detailed  interaction  of  waves  and  turbulence  along  the 
inversion  and  possibly  lead  to  improved  parameterization  of  this  interaction 
in  the  one-dimensional  models. 
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Model  development  during  the  past  year  has  concentrated  on  deriving  a 
hybrid  integral-differential  description  of  the  planetary  boundary  layer  which 
would  allow  approximate  solutions  to  be  obtained  using  far  fewer  numerical 
calculations,  and  on  deriving  a  capability  for  incorporating  precipitation  as 
a  possibility  in  our  two-phase  representation  of  atmospheric  water.  Neither 
of  these  developments  have  reached  the  point  of  being  fully  integrated  into 
our  general  model.  Sections  3  and  4  detail  the  derivation  and  current  status 
of  each. 
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2.  FOG  EVOLUTION  STUDY 


Our  purpose  in  this  section  is  to  demonstrate  the  mechanism  by  which  a 
fog  generated  by  a  decrease  in  surface  temperature  can  relatively  rapidly 
invert  its  internal  temperature  gradient  so  that  the  surface  becomes  warmer 
than  the  air  in  the  fog.  This  mechanism  is,  of  course,  radiative  cooling  of 
the  fog  top  which  then  quickly  moves  down  through  the  fog  layer  by  convective 
overturning.  We  present  a  case  study  of  a  fog  produced  by  warm  airflow  over  a 
cold  surface  illustrating  the  phenomenon;  we  have  considered  the  horizontally 
homogeneous,  time-dependent  problem,  but  this  is  closely  related  to 
steady-state  advection  with  time  corresponding  to  downstream  distance. 

In  order  to  generate  a  cold  water  fog,  we  take  an  initially  fog-free 
boundary  layer,  and  suddenly  reduce  the  surface  temperature.  The  magnitude  of 
the  temperature  change,  AT,  which  is  needed  to  produce  a  fog  depends  on  the 
initial  relative  humidity  of  the  boundary  layer  near  the  surface.  In  the 
absence  of  radiation  effects,  the  surface  layer  analysis  (Reference  10)  can  be 
used  to  give  the  relationship  between  the  critical  AT  and  the  critical 
relative  humidity.  The  relationship  depends  on  the  absolute  temperature  of 
the  surface,  and  the  result  is  shown  graphically  in  Figure  2.1.  It  is  clear 
that  relative  humidities  below  95%  will  require  a  substantial  change  in 
temperature  to  produce  a  fog.  It  should  be  noted  that  radiative  effects  will 
change  this  result  to  an  extent  which  depends  on  the  relative  magnitudes  of 
radiative  and  turbulent  heat  transfer. 

Since  we  are  dealing  with  surface  temperature  changes  of  only  a  few 
degrees,  this  restriction  on  initial  relative  humidity  causes  problems  if  we 
are  trying  to  study  the  fog  evolution  in  isolation.  We  need  to  set  the 
initial  humidity  very  close  to  saturation,  which  means  that  a  fog  is  about  to 
form  even  in  the  absence  of  our  applied  temperature  change  at  the  surface. 
Thus,  the  generation  of  the  fog  depends  as  much  on  the  precise  initial  state 
of  the  boundary  layer  as  on  the  externally  applied  forcing.  We  have  therefore 
not  attempted  any  extensive  study  of  the  fog  evolution  from  different  initial 
states  with  different  surface  temperature  changes,  but  instead  present  a 
single  case  study  and  use  the  integration  as  an  illustrative  example. 

Our  initial  boundary  layer  for  this  case  study  has  93%  relative  humidity 
and  is  only  about  50  m  thick.  The  boundary  layer  was  obtained  by  integrating 
in  time  with  a  constant  geostrophic  wind  equal  to  5  m/sec  and  a  constant  sea 
surface  temperature  for  a  few  days  to  allow  the  humidity  to  increase.  It 
proved  necessary  to  impose  a  significant  subsidence  velocity,  3  cm/sec  at 
1  km,  to  prevent  the  formation  of  a  cloud  layer;  this  is  the  reason  for  our 
relatively  thin  boundary  layer. 

Figure  2.1  shows  that  a  temperature  drop  of  at  least  6°C  is  necessary  to 
generate  a  purely  advective  fog  in  our  initial  boundary  layer  with  its  93% 
relative  humidity.  In  fact,  we  used  AT  =  5°C,  which  means  that  we  rely  on 
radiative  cooling  to  assist  in  the  fog  production,  so  that  the  fog  does  not 
appear  immediately  after  the  surface  temperature  drop.  The  evolution 
following  the  surface  temperature  drop  was  carried  out  under  nocturnal 
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Figure  2.1.  Schematic  showing  the  air  relative  humidity  and 

air-sea  temperature  differences  required  for  pure 
advective  fogs. 
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conditions. 


Figures  2.1  through  2.5  show  the  profiles  of  temperature,  heat  flux, 
liquid  water  content,  and  radiative  cooling  at  several  times  after  the 
change  in  surface  conditions.  We  see  from  Figures  2.1  and  3.1  that  the 
initial  boundary  layer  has  a  warm  surface  and  a  positive  heat  flux  at 
the  bottom;  this  is  required  to  balance  the  radiative  cooling  which  is 
evident  at  t  =  0  in  Figure  2.5 

After  20  minutes,  the  surface  temperature  shows  the  5°C  decrease  from 
290°K  to  285°K,  and  that  the  heat  flux  has  changed  sign  in  the  lowest  15  m. 
There  is  some  cooling  of  the  main  part  of  the  boundary  layer,  by  about  0.5°C, 
due  to  the  radiative  cooling,  and  this  occurs  in  the  absence  of  a  change  in 
surface  temperature.  However,  in  the  latter  case,  the  temperature  only  drops 
by  about  1°C  over  several  hours  and  no  fog  is  formed. 

Figure  2.5  shows  some  radiative  heating  at  the  surface  after  20  minutes; 
but  this  only  extends  about  2  m  vertically,  so  that  the  bulk  of  the  boundary 
layer  is  being  cooled  by  both  turbulent  transfer  and  by  radiation.  This 
causes  the  temperature  to  drop  more  rapidly  than  the  humidity  (which  is  only 
reduced  by  turbulent  transfer)  so  that  the  air  eventually  saturates  and  a  fog 
forms  around  t  =  45  mins.  At  this  stage  the  air  has  cooled,  so  that  the 
surface  heat  flux  is  reduced  in  magnitude  from  its  value  at  20  minutes,  but 
the  air  is  still  warmer  than  the  surface.  The  fog  extends  about  12  m 
vertically  at  this  stage  (see  Figure  2.4).  Also  at  t  =  45  min,  Figure  2.5 
shows  the  beginning  of  the  increased  radiation  from  the  top  of  the  fog. 

After  t  =  45  min  the  fog  develops  in  depth  and  intensity,  and  the  air 
temperatures  continue  to  fall  due  to  radiative  cooling  from  the  fog  top.  The 
heat  flux  profiles  show  that  the  lower  part  of  the  fog  is  cooled  by  turbulent 
transfer.  Shortly  after  t  =  80  m,  the  air  temperatures  drop  below  the  surface 
temperature,  and  we  have  a  warm  surface  fog  thereafter.  The  development 
continues  after  this  time  with  the  depth  of  the  fog  layer  increasing  and  both 
turbulent  and  radiative  heat  fluxes  also  increasing. 

Thus,  in  this  particular  case,  a  cold  water  fog  develops  roughly  40 
minutes  after  the  change  in  surface  temperature,  and  persists  as  a  cold-water 
fog  for  a  further  40  minutes;  after  this  time  it  is  converted  to  a  warm-water 
fog  and  continues  to  deepen  and  intensify.  One  may  expect  that  fogs  initiated 
by  smaller  drops  in  surface  temperature  will  transition  to  a  warm  water  fog 
more  quickly.  For  a  cold-water  fog  to  persist  it  appears  necessary  to  have 
both  relatively  strong  winds  and  a  drop  in  surface  temperature  which  is 
stronger  than  that  indicated  in  Figure  2.1  for  the  particular  ambient 
relative  humidity. 

Figure  2.4  shows  that  the  liquid  water  content  in  the  calculated  fog 
increases  rapidly  up  to  an  equilibrium  level  of  about  lg/kg.  This  value 
appears  to  be  significantly  higher  than  atmospheric  measurements  which 
typically  give  liquid  water  contents  of  0.1  -  0.2  g/kg.  It  seems  likely  that 
the  reason  for  this  discrepancy  is  the  absence  of  any  mechanism  for  removal  of 
liquid  water  through  gravitational  settling  in  the  model.  In  measurements 
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Figure  2.2.  Evolution  of  air  virtual  potential  temperature 
following  a  drop  of  5°C  in  the  sea  surface 
tempera ture. 
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45  minutes 
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Figure  2.3.  Evolution  of  the  turbulent  heat  flux  following  a  drop  of  5  C  in  the  sea  surface 
tempera ture. 


Figure  2.4.  Evolution  of  the  liquid  water  content  following  a  drop 
of  5°C  in  the  sea  surface  temperature. 
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Figure  2.5  Evolution  of  the  radiative  cooling  rate  following  a 
drop  of  5°C  in  the  sea  surface  temperature. 


17 


over  land,  (Reference  19)  concluded  that  only  a  small  fraction  of  the  liquid 
water  which  condensed  actually  remained  in  the  boundary  layer.  The  balance  was 
presumed  to  have  been  deposited  on  the  ground  by  gravitational  setting.  In 
our  model,  we  assume  a  fixed  droplet  size  spectrum  for  the  purposes  of  the 
radiation  calculation,  and  ignore  gravitational  settling  since  this  is 
negligible  for  our  assumed  spectrum.  A  significant  improvement  would  probably 
be  obtained  from  a  cloud  physics'  model  which  accounted  for  droplet  growth 
within  the  fog  and  allowed  the  heavier  drops  to  fall  out.  These  processes 
have  not  been  unambiguously  identified  as  dominant  controlling  mechanisms  on 
the  liquid  water  content  in  atmospheric  fog  studies,  but  they  are  the  most 
1 ikely  candidates. 
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3.  HYBRID  INTEGRAL-DIFFERENTIAL  DESCRIPTION  OF  THE 
PLANETARY  BOUNDARY  LAYER 
USING  SECOND-ORDER  CLOSURE  TURBULENCE  THEORY 


3.1.  Introduction 

In  this  work  we  describe  a  hybrid  integral-differential  procedure  for  the 
prediction  of  the  dynamics  of  the  horizontally  homogeneous  planetary  boundary 
layer  (PBL)  according  to  a  second-order  closure  theory  of  turbulence.  This 
procedure  is  aimed  at  enhancing  the  quality  of  the  usual  finite  difference 
solution  of  these  equations  by  the  incorporation  of  integral  constraints.  In 
particular,  the  intent  is  that  a  finite  difference  solution  utilizing  on  the 
order  of  five  grid  cells  should  provide  adequate  accuracy  with  the  hybrid 
method.  The  result  would  be  a  computationally  efficient  procedure  which  would 
be  useful  in  operational  applications.  Such  a  method  used  as  the  basis  for 
the  inhomogeneous  (three-dimensional)  PBL  would  offer  similar  computational 
advantage  and  efficiency. 

In  Section  3.2  we  present  the  basic  concept  of  the  method.  In 
Section  3.3  we  indicate  the  manner  in  which  integral  constraints  are 
incorporated  with  the  finite-difference  solution.  In  Section  3.4,  we 
illustrate  the  method. 


3.2.  Concept  of  the  Hybrid  Method 


The  second-order  closure  theory  of  the  PBL  developed  and  in  use  at 
A.R.A.P.  (References  2,  7,  10)  has  proved  quite  successful  in  describing  the 
dynamics  of  the  PBL  under  most  general  steady  state  and  transient  situations 
and  from  stable  to  unstable  conditions,  including  the  presence  of  a  capping 
inversion  layer.  Even  for  the  horizontally  homogeneous  layer,  the 
computational  volume  for  executing  these  descriptions  is  still  considerable, 
requiring  30  to  50  vertical  grid  levels  for  adequate  resolution,  particularly 
if  an  inversion  is  present.  Many  of  the  principal  quantities  of  interest  in 
the  PBL  are  global,  or  surface  quantities,  rather  than  detailed  local  interior 
quantities.  These  key  quantities  include  the  surface  fluxes  of  momentum, 
heat,  and  species,  as  well  as  the  total  boundary. 1 ayer  depth  and  the  cross 
isobaric  wind  angle.  As  such  they  are  described  by  the  integral  forms  of  the 
equations  of  motion;  however,  these  integral  forms  contain  integrals  over  the 
profiles  of  the  mean  and  turbulent  field  variables. 

There  are  two  regions  of  sharp  gradients  in  the  general  PBL:  the  surface 
layer  and  the  inversion  layer.  It  is  therefore  useful  to  treat  these  regions 
as  "integral"  regions  in  which  analytical  or  approximate  forms  of  the  profiles 
are  utilized  in  integral  forms  of  the  equations  of  motion.  The  "outer  layer" 
which  exists  between  the  surface  layer  and  inversion  layer  may  then  be  treated 
either  in  integral  ("single  layer")  fashion  or  with  a  full  finite  difference 
treatment. 
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The  hybrid  integral -differential  method  combines  integral  constraint 
equations  with  the  full  finite  difference  equations  which  describe  the  outer 
layer.  The  method  can  thus  be  used  either  with  a  large  number  of  finite 
difference  levels  or  with  a  sparse  number  of  levels  with  a  "good"  solution  for 
the  boundary  layer  parameters  still  achievable  in  the  limit  in  which  the 
finite  difference  levels  consist  of  no  points,  i.e.,  the  description  is  purely 
integral.  In  this  limit,  the  method  bears  some  resemblance  to  various  "layer" 
parameterization  procedures  (References  20,  21,  22).  However,  there  is  an 
essential  difference  between  these  layer  parameterization  schemes  and  the 
hybrid  procedure  we  describe  here.  Mathematical  approximation  for  simplicity 
of  solution  and  physical  modelling  interplay  together  in  the  formulation  of 
the  layer  parameterization  models  cited  above.  The  hybrid  procedure  described 
here  including  its  integral  constraints  is  rigorously  based  upon  the  second 
order  closure  theory  of  turbulence.  Hence,  turbulence  modelling  issues  are 
confined  to  the  validity  of  the  general  second  order  closure  theory.  The 
hybrid  procedure  is  directed  to  the  representation  of  this  system  of  equations 
and  the  method  of  solution  of  the  system.  Thus,  in  the  hybrid  method 
described  here,  solution  approximation  issues  are  separated  from  turbulence 
modelling  issues.  The  capability  for  continuous  transition  from  a  purely 
integral  description  to  a  fully  differential  description  is  one  of  the 
features  we  have  attempted  to  incorporate  in  the  hybrid  procedure.  Of  major 
interest  is  the  intent  that  the  integral  constraints  which  are  part  of  the 
hybrid  method  will  significantly  improve  the  solution  of  the  PBL  for  cases  in 
which  a  sparse  grid  of  only  four  to  six  finite  difference  levels  is  used  over 
that  which  would  be  obtained  in  absence  of  these  integral  constraints. 

A  full  second-order  closure  PBL  model  which  could  adequately  perform  with 
only  four  to  six  finite  difference  cells  would  be  computationally  efficient 
and  would  possess  advantages  for  operational  implementation  in  various 
applications  as  well  as  providing  an  economical  basis  for  fully 
three-dimensional  (horizontally  inhomogeneous)  PBL  descriptions. 

The  full  range  of  PBL  behavior  ranging  from  stable  to  unstable 
including  those  capped  by  an  in  inversion  layer  may,  in  principle,  be 
treated  with  the  hybrid  method.  Large  scale  divergence,  humidity 
effects,  and  radiative  transport  represent  further  key  processes 
which  require  inclusion  in  the  hybrid  method.  We  have  selected  the 
following  reasonable  steps  of  development  and  evaluation  of  the  hybrid 
method  for  the  I-D  homogeneous  PBL. 

1 .  Neutral  PBL 

2.  General  PBL  (stratified) 

3.  General  PBL  including  inversion  layer  and  large  scale 

divergence 

4.  General  Moist  PBL  (including  humidity  and  radiative 

transport) 

In  the  present  report  steps  1  and  2  have  been  completed  and  sample  results 
are  presented.  These  results  are  worthy  enough  to  encourage  us  to  begin 
steps  3,  4.  The  major  activities  in  these  steps  involve  the  development 
of  the  integral  equations  for  the  inversion  layer  as  well  as  integral 
equations  for  the  prediction  of  the  humidity  boundary  layer  thickness 
and  fluxes  and  the  cloud  base  and  cloud  top  elevations. 
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3.3  Integral  Constraints  and  Their  Incorporation  Into  the 
Hybrid  Method 


3.3.1  Governing  Equations  of  the  Horizontally  Homogeneous  PBL 

We  consider  a  horizontally  homogeneous  PBL  in  coordinates  x,y,z  aligned 
with  z  normal  to  the  surface  and  x,y  as  coordinates  in  the  plane  perpendicular 
to  z.  Let  0  =  (J(U,V)  be  the  velocity  field  in  the  plane  and  let  (UQ,Va)  be 
the  geostrophic  wind  velocity  components  defining  the  pressure  gradients.  The 
momentum  equations  in  the  (x,y)  plane  may  be  expressed  as 


gt  -  U)  =  ”  (UW)  -  f  (V  -  Vg)  +  U0 


(3.1) 


at 


(v.  -  V)  =  —  (vw)  +  f  (u  -  Uq)  +  V0 

3  Z  ^ 


(3.2) 


In  the  above,  f  is  the  Coriolis  parameter,  uw,  vw  are  the  stress  components 
and  -  OJU^VJ  is  the  velocity  vector  in  the  inviscid  region  above  the 

boundary  layer.  The  corresponding  mean  thermal  energy  equation  for  the 
virtual  potential  temperatuure  o  is 


_JL  v 

(Wb) 

az 


Q  +  Q 


(3.3) 


where  Wb  is  the  vertical  turbulent  heat  flux  and  Q  is  the  thermal  energy 
source  term.  The  inviscid  region  forms  of  Eqs.  (3.1)  through  (3.3)  are 


au« 

at 


f(V„-  Vgj  -  0O 


(3.4) 


aVco 

at 


+  f  (U.  -  u  ) 

y  ao 


(3.5) 
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=  Qco  (3.6) 

where  Uq^,  Vgoo  define  the  pressure  gradients  outside  the  boundary  layer.  In 
what  follows  we  shall  assume  Ug  =  U9oo,  Vg  =  Vg„.  Thus  Eqs.  (3.1)  through 
(3.3)  may  be  expressed  as 


n  (u~  -  u) 

i  (v~  -  v> 

7t  (6>  -  0) 


9  (uw)  +  f  (V.  -  V) 

o  Z 

~  (vw)  -  f  (U.  -  U) 

7“  (wa)  -  (Qoo  +  Q) 

o  2 


(3.7) 

(3.8) 

(3.9) 


The  turbulent  moment  equations  at  second-order  closure  level  which  we  shall 
employ  include  those  for  the  Reynolds  stresses,  heat  flux,  temperature 
variance,  and  the  turbulent  scale  equation.  We  do  not  repeat  the  full  set 
of  turbulence  equations  here,  but  refer  the  reader  to  Reference  7.  The 
two  equations  from  the  turbulence  set  which  we  will  repeat  here  because  of 
their  use  in  the  integral  developments  are  the  turbulent  kinetic  energy 
equation  for  the  turbulent  kinetic  energy  0 /2)q2  and  the  turbulent  scale 
equation  for  the  scale  A.  These  are 


a 

at 


—  all  —  aV  g  — 
-  uw  -  vw  —  +  wo 
3z  dz  0O 


(3.10) 
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(3.11) 


9A 

aT 


( —  3A  —  3A 

luw - +  vw - 

\  9Z  3Z 


) 


+  0.075  q  +  0.3  — - 
M  dz 


-MU  (iai\  +  0.84JS 

q  \  9z/  qz  t0 


Eqs.  (3.7)  through  (3.10)  along  with  the  remaining  second-order  closure 
equations  form  the  basis  for  the  subsequent  development. 

The  mean  defect  kinetic  energy  equation  from  Eqs.  (3.7)  and  (3.8J  is 


3£  (m) 

8  t 


(Uco  -  u)  —  (uw)  +  (V.  -  V)  —  (vw) 
<*z  8Z 


where 


E(m) 


1 

2 


U 


2 


The  total  energy  equation  derived  from  Eq.  (3.10)  and  the  above  equation 


3E 

at 


a 

9Z 


uw  (Uoo  -  U)  +  vw  ( Voo  -  V ) 


°-3  tz  [qA  tz  (i q")] 


•  12) 
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3.3.2  Integral  Parameters  and  Integral  Constraints 


The  important  integral  (global  as  opposed  to  local)  parameters  of  the 
neutral  PBL  are  readily  cited  as  the  surface  shear  stress  or  shearing  velocity 
components  (u*,v*)  and  the  overall  boundary  layer  depth  (6).  In  addition  to 
these,  integral  parameters  may  be  generated  by  simply  taking  moments  over  the 
mean  and  turbulent  fluid  field  variables.  The  first  moments  are  the  average 
velocity  components: 


6 

<u>  =  6o1  J  dz 

zo 


<V> 


0 

*7 


V(z)  dz 


In  the  above  equations  we  have  6q  =  (6  -  zo)  where  Zo  is  the  effective 
roughness  height.  These  average  velocities  are  equally  expressible  in  terms 
of  the  displacement  thickness  6u>c>v  defined  as 


6 u  =  lUj"1  J  (Uco-U)dz  6V  =  |Ujy  (V.-tf)  dz  (3.13) 

Z0  z0 

where  U  =  U(U,V)  is  the  local  velocity  vector,  is  the  velocity  vector  in 
the  inviscid  region  above  the  boundary  layer  and  (U^.V^)  are  the  velocity 
components  in  this  inviscid  region. 

In  a  stable  PBL,  the  surface  heat  flux  expressed  in  terms  of  ti*  will 
enter  the  set  of  integral  parameters  as  well  as  the  average  temperature  <e>, 
or  equivalently,  a  thermal  energy  thickness  6d  defined  as 


6e  =  e;‘  /  (e-  - e)  dz 

z0 


(3.14) 


here  0  is  a  reference  temperature  and  6^  is  the  virtual  potential  tempera¬ 
ture  iFI  the  inviscid  zone  at  the  boundary  layer  edge.  If  other  energy 
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source  terms  appear  in  the  thermal  energy  equation  (such  as  radiative  source 
terms)  the  temperature  requires  more  detailed  definition  appropriate  to  the 
particular  source  term  in  question.  Note  that  6^  may  be  a  function  of  Z. 

In  an  unstable  PBL,  the  above  parameters  also  enter;  however  in  addition 
a  significant  set  of  integral  (over  the  inversion  layer)  parameters  describing 
the  inversion  layer  will  also  enter.  We  do  not  take  up  the  additional 
integral  parameters  and  constraint  equations  for  the  inversion  layer  in  this 
report  as  described  on  page  20. 

The  integral  constraint  equations  corresponding  to  the  foregoing  integral 
parameters  are  the  first  spatial  moments  of  Eqs.  (3.7)  through  (3.9).  These 
are  (assuming  zq  independent  of  t): 


£(|9J«„>  *  +  f6v|0„[ 


~  (|u„|6v)  =  vj  -  f6u|U„| 


(er60>  =  *  se 


(3.15) 

(3.16) 

(3.17) 


The  surface  shgar  velocity  vector  is  u* 
v*  =  -(vw)0,  |u*|«*  =  -(wb)  and  (  )0 
energy  source  term  se  is  defined  by 


=  u|  (u*,v|)  where  u^  = 
denotes  a  surface  value 


-(uw)0. 

The  thermal 


s  0 


J  [Q®  -  Q(z )]  dz 


The  next  moments  involve  products  of  the  velocity.  The  defect  mean 
kinetic  energy  thickness  6^m)  is  defined  as 


6 


Zo 


i 


2 


dz 


(3.18) 
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where  Er  is  a  reference  energy  which  requires  specification  to  complete 
definition.  The  turbulent  kinetic  energy  thickness  may  be  defined  as 


6  (t) 
E 


q2  dZ 


where  q  is  the  RMS  turbulence  velocity.  One  can  also  consider  the  total 
energy  thickness  <$£  defined  as 


6 


E 


(3.19) 


Since  it  is  the  turbulent  viscous  field  which  underlies  the  defect  mean 
kinetic  energy  as  well  as  the  turbulent  kinetic  energy,  it  is  appropriate  to 
select  the  reference  energy  Er  as  the  surface  turbulent  kinetic  energy 
expressed  in  terms  of  the  surface  RMS  turbulence  velocity  q0  for  neutral  or 
stable  boundary  layers: 


E 


r 


(3.30) 


We  describe  the  method  of  selecting  Er  in  the  case  of  unstable  boundary  layers 
in  Section  3.6.  The  integral  total  energy  equation  has  the  form 


at 


r\  a 

-  (E  6C)  =  u  Uo  +  vlt 
r  E  *  ★ 
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wo  dz 


1 
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6 

[  3.  Iq2  dz  (3.31) 

J  A  L 
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The  characteristic  average  scale  in  the  boundary  layer  may  be  defined  as 


(3.32) 
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The  integral  equation  governing  z  is  then  the  integrated  form  of  the  scale 
equation  (3.11): 


—  (z6  )  =  0.35 
9t  0 


/ 


—  (uw  —  +  vw  — )  dz  +  0.075 
q2  az  az 


6 

/qdz 

20 


-  0.375  f6  I  (-1  qA\2  dz  -  0.8  _£  wu  dz  (3.23) 

/  q  la  z  /  J  q2  T0 

It  can  be  seen  that  the  Eqs.  (3.15  -  3.17)  and  (3.21  -  3.23)  form  a 

closed  system  for  the  integral  parameters  <5ui  6V,  6e,  6^,  z  provided  the 
fluxes  u*,  v*,  0*  and  the  integrals  appearing  in  Eqs.  (3.21)  and  (3.23)  can 
be  expressed  in  terms  of  these  quantities. 
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3.3.3  Structure  Function  Approach  to  the  Surface  Flux  Laws  and  the 
Integral  Parameters 


We  now  take  up  the  method  of  closure  of  the  hybrid  procedure  which 
unifies  the  integral  constraints  with  the  fully  differential  procedure.  The 
unifying  vehicle  is  a  system  of  structure  functions.  We  consider  first  the 
representation  of  the  surface  fluxes  in  terms  of  the  integral  parameters. 

In  the  absence  of  stratification  effects,  the  velocities  and  turbulent 
momentum  fluxes  in  the  surface  layer  in  a  right  handed  orthogonal  coordinate 
system  with  the  U  component  aligned  parallel  to  the  outer  flow  velocity 
may  be  represented  as 

uw  =  -  | lT* j  cose  (3.24) 


vw 


sine  (1  -cj>- ) 


(3.25) 


U 


K 


In 


(z/zq)cos3 


(3.26) 


V  = 


(z/2q )si n3(l  -4>2) 


(3.27) 


where  u*  =  (u£,v*)  is  the  vector  of  surface  stress  components  lying  in  the 
plane  of  motion.  Here  uw  ,  vw  are  the  turbulent  momentum  fluxes  parallel 
and  perpendicular  to  the  outer  flow  velocity  while  U,  V  are  the  ^ 
corresponding  velocity  components.  The  angle  between  the  local  velocity  U 
and  the  outer  flow  is  denoted  as  6  .  The  absolute  magnitude  of  the 

surface  shear  stress  components  lying  in  the  plane  of  the  flow  is  |u£|  . 

The  functions  <f>](z)  ,  ^(z)  are  representations  of  the  effects  of  the 
Coriolis  force-induced  pressure  gradient  in  forcing  vw  and  V  to  deviate 
from  constant  and  logarithmic  values  respectively  which  is  the  first 
manifestation  of  the  pressure  gradient  effect  on  the  surface  layer  for 
z  >  zQ  ,  (Reference  23): 

f|ffj  |u*| 

(z)  =  r  ? -  [z  -  zn - T —  cose  zy(z)]  (3.28) 

1  iu*l2sine  0  K|uj 
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fin 
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u*|Sln6  ln(z/z  )  1  0  0 


[■ 


—7 — cose  zy(z)  -  [z  ln(z/z  )  -  (z 
k|U  |  1 


-z°)jD 


(3.29) 


with 


y(z)  =  1n(z/zo)  -  (1-zQ/z) 


Now  let  am  denote  the  angle  of  the  outer  flow  velocity  with  respect 
to  the  x  axis  of  an  arbitrary  coordinate  system  fixed  to  the  earth.  Let  a 
denote  the  angle  of  the  local  velocity  ti  makes  with  the  same  coordinate 
system.  The  angle  3  is  related  to  ,  a  as 


3  =  a 


a 


00 


The  surface  layer  momentum  flux  and  velocity  components  in  this  general 
coordinate  system  in  which  we  shall  formulate  the  procedure  are  given  by 


uw  =  - 


[cosa  +  sin3  sina^  <f>-|(z)] 


(3.30) 


vw 


2 


[sina  -  sin3  cosa^  $-|(z)] 


(3.31) 


U 


l'n(z/z0) 


[cosa  +  sin3  sina^  (^(z)] 


(3.32) 


V 


ln(z/zQ) 


[sina  -  sin3  cosa^  ^(z)] 


(3.33) 


Consider  z  =  6S  as  the  height  of  the  surface  layer.  Then  at  this 
level  we  have  from  Eqs.  (3.32)  and  (3.33) 


<\\\ 

$(6s)ln(6s/ZQ) 


(3.34) 
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where 


$(6S)  =  {1  +  sin23  <}>2(6s)[(J)2(6s)  -  2]}1/2 


Let  us  now  compare  this  result  with  the  Monin-Obukhov  surface  layer  similarity 
theory  for  stratified  flows  which  (neglecting  Coriolis  effects)  has  the  form 


ln(6s/zo)  +  ia 


u 


(3.35) 


where  yu  =  y,j(Ss/L)  is  a  Monin-Obukhov  similarity  function  for  momentum 
transfer  and  has  different  forms  depending  upon  whether  the  Monin-Obukhov 
length  L  is  greater  or  less  than  zero.  The  Monin-Obukhov  length  is  defined 
as 


T  I 


->-2  i 


L  .-2- 


gK 


We  observe  that  the  form  Eq.  (3.34)  which  includes  first  order  corrections 
for  the  effects  of  the  Coriolis  forces  in  the  surface  layer  was  derived  with 
the  condition  that  u*  and  the  RMS  turbulence  velocity  qq  were  approximately 
constant  in  the  surface  layer.  On  the  other  hand,  the  Mjgnm-Obukhov 
similarity  form  (3.35)  is  based  on  the  condition  that  |u*|  is  constant  in 
the  surface  layer.  For  flows  in  which  the  angle  3  is  not  too  large,  it  may 
then  be  possible  to  obtain  a  general  "extended"  surface  layer  resulting  in 
the  form 


where 


D 


u 


H  *(«s/z0)ln(Ss/20)  +  ug 


(3.36) 


(3.37) 


Although  we  have  not  proved  it  here,  we  conjecture  that  the  form  Eq.  (3.37) 
is  correct  to  first  order  in  the  surface  layer  expansion  parameter 


VLc 


fit)  |  6 

1  co  I  $ 

|u*l |sin3| 


(3.38) 
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where  Eq.  (3.38)  defines  the  Coriolis  surface  length  Lc  which  can  be 
interpreted  as  the  height  at_^which  Coriolis  force  induced  pressure  gradients 
disturb  the  uniformity  of  |u*|  and  the  pure  logarithmic  form  of  the 
velocity  variation  in  the  surface  layer. 

We  note  that  by  choosing  6S  «  |L|  ,  6S  «  Lc  we  may  always  render 
the  form  of  Du  in  Eq.  (3.37)  in  the  pure  logarithmic  form.  There  is  a 
usefulness,  however,  in  allowing  the  surface  layer  to  be  as  thick  as  possible 
with  consideration  for  the  validity  of  Eqs.  (3.30)  -  (3.33)  in  application 
of  the  hybrid  procedure  with  a  space-grid  finite  difference  procedure. 
Incorporation  of  virtually  the  full  logarithmic  layer  below  z  =  6S  then 
allows  the  finite  difference  procedure  a  better  resolution  of  the  more  linear 
region  for  z  6S  .  In  the  absence  of  stratification,  Eqs.  (3.30)  -  (3.33) 
are  quite  accurate  even  for  5S  >  Lc  yielding  results  that  are  within  10%  of 
the  exact  solution  for  the  steady  state  PBL  at  <5S  =  1 0  Lc  . 

The  counterparts  to  Eqs.  (3.30)  -  (3.33)  for  heat  flux  and  temperature 
distribution  in  the  surface  layer  are 


w0  =  - 


(3.39) 


0  "  0o  "  0*PRDe/|< 
from  which  we  may  write 


0*  - 


PR°0 


(e,  -  eJ 


(3.40) 


(3.41) 


where  0  is  the  virtual  potential  temperature  at  z  =  6S  ,  0Q  is  the 
surface  temperature,  0r  is  a  reference  temperature,  and  Pp  is  the 
turbulent  Prandtl  number.  The  function  Dg  like  Dy  consists  of  a 
logarithmic  portion  and  a  stratification  portion  embedded  in  the  Monin-Obukhov 
similarity  function  p0  : 


De  -  i"(yzo>  ♦  WL> 


(3.42) 


To  relate  the  fluxes  |u*|,  0*  ,  and  the 

to  the  integral  parameters  <5U  ,  6V  ,  6g 
0S-0O  in  terms  of  the  structure  functions 


surface  layer  angle  a  (or  8) 
we  express  ,UC  ,  Yc  >  and 
stuT  ,  s(V),Sste)S  as  follows 
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0 

(3.43) 

0 

(3.44) 

=  S^<0-0  > 

O  0  0 

(3.45) 

where  the  (  )'  indicates  representation  in  the  coordinate  system  aligned  with 
ft  and  the  <  >0  indicate  a  spatial  average  over  the  outer  part  of  the 
boundary  layer: 


<(  )>0  -  («-«,)■ 


/ 


(  )dz 


Tlpe^e  av^r^ges  may  in  turn  be  related  to  the  displacement  thicknesses  6^  , 

as 


$57¥*FB“. 

Lio  -  c°saoo  6io)/ 

*1  I  /oinn, 


-  sina  ) 

SO  oo  v  so 

(3.46) 

cosa  6  f,0  ^  /  6  ) 
oo  v  so 

(3.47) 

0°>/6so) 

(3.48) 

i  the  structure  functions 

,  sty) .  s 

parameters 


ie  Tiux  i avii  wiiiLii.reiai 

6 (° )  ,  ,  6a0'  as 


K|ftq 


(3.49) 


(3.50) 
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(3.51  ) 


e  =  a  -  a  =  tan-1  (V'/l)' ) 

00  S  S 


with  U'  ,  V'  and  0  -0  given  as 
s  s  °°  s  3 


Ul  =  It)  |S(U)(1  -  cosa  6(o)/6c  -  sina  6^/6  ) 

S  1  ool  oo  u  SO  OO  v  so 


V'  =  |tMs^(sina  -  cosa  S^^/S) 

s  *  ool  00  u  'so  oo  v  SO 


|tis|  =  /(up2  +  (vp2 


(0  -0  ) 
v  OO  $  1 


[  0~‘9O 
\  er 


The  "outer"  thicknesses  6^°^  ,  6^°^  ,  6q0^  are  related  to  the 
thicknesses  6  ,  6„  ,  6Q  as 

U  V  u 


where  the  inner  thicknesses  6^  are  related  to  |u*|  ,3,0, 

J 


=  (6  /s  )|ti  r1 

u  SO  0  1  00 1 


(u  -U)dz 


i*1*  =  <«„/«„) i®  r1 

V  so  o  '  °°l 


(V  -v) dz 


(3.52) 

(3.53) 

(3.54) 

(3.55) 

total 

through 

(3.56) 

(3.57) 
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(eco~0)dz 


(3.58) 


(i)  _ 


(6  /6  )0 
v  so7  o'  r 


s 

/ 


where  U  ,  V  ,  0  in  Eqs .  (3.56)  -  (3.58)  are  given  in  terms  of  the  surface 
layer  expressions  (3.32),  (3.33),  and  (3.40). 

Equations  (3.49)  -  (3.55)  are  the  key  results  of  this  section.  They 
provide  the  surface  flux  laws  which  relate  u*  and  9*  to  the  integral 
parameters  6U  ,  6V  .  Sq  and  in  so  doing  naturally  introduce  the  structure 
functions  S(U),  S^K  s(0)  which  are  obtained  from  detailed  solution  of 
the  differential  equations.  Equations  (3.36)  and  (3.41)  are  exact  statements 
provided  the  structure  functions  s(U)  ,  s(V)  ,  s(9)  are  known.  These 
equations  relate  the  surface  fluxes  u*  ,  v*  ,  0*  to  the  integral  thickness 
<5U  ,  Sv  ,  6q  .  As  such,  they  may  be  considered  friction  and  heat  flux  laws 
for  integral  PBL  description.  We  now  observe  that  if  the  detailed  profiles 
U(z),  V(z),  0(z)  are  known,  the  structure  functions  ,  s(V)  ,  s(0' 

may  be  directly  calculated  (with  specification  of  the  thickness  6S).  Hence, 
these  structure  functions  may  be  calculated  in  terms  of  the  profiles  generated 
by  the  finite  difference  solution.  It  may  be  further  observed  that  the 
integral  constraints  furnish  u*  ,  v*  ,  0*  in  terms  of  the  integral  thickness 
6U  ,  Sy  ,  6q  .  Hence,  these  surfaces  fluxes  together  with  the  surface  layer 
functions  Eqs.  (3.52,  3.33,  and  3.40),  determine  the  boundary  data  at  z  =  6S 
for  the  finite  difference  description  of  the  domain  6S  <  z  <  6  .  As  such, 
the  derivative  boundary  condition  normally  required  at  the  top  of  the  surface 
layer  z  =  6S  for  the  finite-difference  equations  is  dispensed  with.  The 
derivative  boundary  conditions  are  thus  replaced  by  the  integral  constraints. 
The  interest  is  that  this  procedure  for  treating  the  surface  conditions  should 
significantly  improve  the  quality  of  the  overall  finite  difference  description 
when  a  sparse  set  of  grid  levels  is  utilized  consisting  of  perhaps  only  four 
to  six  grid  cells. 
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3.3.4  Transfer  Coefficients 


We  may  define  transfer  coefficients  Cf  ,  Cq  for  momentum  and  heat 
transfer  to  the  surface  as 


I  +2 1 

u* 


1112 


(3.59) 


ce  = 


I 


u 

OO  I 


(3.60) 


The  momentum  and  heat  fluxes  are  then  expressed  as 


J*  =  cf 


2 

'*  =  cf 


\t  1 

2 

cosa 

(3.61 ) 

1  00  • 

It)  1 

2 

1  sina 

(3.62) 

1  00  1 

C0 

8=+ 

CD 

(3.63) 

It  is  also  useful  to  define  differential  momentum  and  heat  transfer  coefficients 
which  characterize  the  rates  of  change  of  u*  ,  v*  ,  |u*|0*  with  respect  to 
the  integral  parameters  5U  ,  6„  ,  60  (with  all  other  quantities  held  fixed). 
Let  6a  represent  the  vector  of  thicknesses  6U  ,  6V  ,  60  and  ca  the  vector 
of  fluxes  CfCOSa  ,  CfSina  ,  c0  : 

l6u\ 

6  =  6..  ] 
a  V 

w 


c  = 

a 


c^cosa 


c^sina 


\ 


(3.64) 


The  differential  exchange  coefficients 


C  0  are  then  defined  as 

CXp 
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(3.65) 


Ca3  =  6o  3V963 


For  the  case  in  which  =  0  these  exchange  coefficients  have  the  form 


'12 


'21 


'22 


rS(U) 


'f  ft 


:(U) 


(V) 


(V) 


1+3A  \ 

cos2B  +  sin2B 

(  1 +3Au  \ 

^2  — -  lj  singcosB 

'  1  +3AV  \ 

^2  — j - 1  I  cos3sin3 

r  / 1  +3i  \  - 

■j —  I  sin  6  +  cos  6 


(3.66) 


'13 


'23 


c  |  — -§ —  \  s(0) 
fie  -0  1  b 
\  s  0 


A 

_ L 

A 


2  -r^- 


cos3 


sin3 


(3.67) 


'31 


'32 


1  -3Ar 


\  =  c 


6  rfr 


S(U)  cos3 


sin3 


(3.68) 


C33  "  C0  10  -6 


s  o 


•  (e)  (  1-Ae 


(3.69) 


In  the  foregoing 


A  .  i- 


u  Du  3L 


A  .JL 

8  °8 


9L 
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1 


4  =  1  +  2iu  -  40 


In  the  stable  limit  where  uu  ~  ,  y0  ~  j-  the  functions  Au  ,  A^  become 


Au  =  -"t 


Ae  d. 


and  obey 


-  1  *  Ay  <  0 


1  <  A0  <  0 


wi  th 


Au  ,  Ae  ->  °  as  L  - 


A  ,  A0  -»■  -1  as  L  -»■  0 


We  may  note  that  as  L 


0  and  the  coefficients 


-ly  •  «  anu  me  V.UCI  i  IV.  ICHW  C]  3  ,  C23  0 

The  momentum  flux  then  becomes  decoupled  from  the  heat  flux.  On  the  other 
hand,  for  finite  L 
case  where  a  =  0  , 
given  by 


the  heat  flux  and  momentum  flux  are  coupled.  In  the 
the  Brunt  Vaisaila  frequency  in  the  surface  layer  is 


O), 


(S) 

BV 


& 


3C31 


ft  |/6 

oo  I  '  r 


(3.70) 


When  L  <  0  , 
frequency  is 


the  product 
imagi nary. 


C13C31 


is  less  than  zero  and  the  Brunt  Vaisaila 
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3.3.5  Surface  Layer  Depth 


Let  us  now  consider  the  surface  layer  depth  6 
determined  by-  the  strength  of  the  terms  in  Eqs.  (3. 
condition 


This  depth  is  rigorously 
3.9)  which  disturb  the 


vw 


W0  ) 


0 


and  the  validity  of  the  surface  layer  solutions  Eqs.  (3.24  -  3.27).  There  are 
three  general  effects:  (1)  the  Coriolis  effect,  (2)  unsteady  effects,  and 
(3)  thermal  source  terms.  From  the  work  of  Reference  23  it  can  be  shown 
that  the  Coriolis  effect  leads  to  6$  determined  by 

6$  <  Lc  =  I sine|  |u*|/(|f|  jtfj)  (3.71) 


Let  us  now  consider  the  unsteady  effects  characterized  by  unsteady  forcing 
frequencies  f^°°)  ,  f\°° '  .  Such  effects  lead  to  conditions  of  the  form 


«s  S  iSfl/tlfJf’l  I'Ll)  «s  <  |u.|  |e*|/(|fg")|er)  (3.72) 

Under  all  conditions  the  surface  layer  depth  6S  should  be  limited  by 

6S  <  | L |  (3.73) 

We  remark  that  although  the  expedient  of  setting 

6$  =  e  60  (3.74) 


(where  e  is  some  small  fraction)  is  simple  and  attractive,  such  a  procedure 
is  not  necessarily  consistent  with  the  definition  of  the  surface  layer  as  that 
region  in  which  the  solutions  Eqs.  (3.24  -  3.27)  are  valid.  Although  Eq. 

(3.74)  would  likely  lead  to  a  determination  of  6S  consistent  with  Eqs.  (3.24  - 
3.27),  the  conditions  (3.71),  (3.72),  and  (3.73)  are  the  more  rigorous  conditions 
for  the  determination  of  . 
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3.3.6  Energy  Thickness  and  the  Dynamic  Equation  for  the  PBL  Thickness 


We  now  consider  the  development  of  the  integral  energy  equation  in  terms 
of  integral  structure  functions.  We  consider  the  choice  of  a  characteristic 
energy  Er  and  characteristic  RMS  turbulence  velocity  qr  .  In  a  stable 
PBL  (L  >  0)  we  select  qr  as  the  surface  value  q0  .  In  the  case  of  an 
unstable  PBL  (L  <  0)  the  surface  turbulence  level  does  not  characterize  the 
average  turbulence  levels  within  the  PBL  because  of  the  buoyant  production 
of  turbulence.  In  this  case,  we  select  qr  as 


Hence,  the  reference  turbulence  and  energy  levels  are  given  by 


L  >  0 


(3.75) 


q 


r 


L  <  0 


E, 


(3.76) 


We  now  represent  the  integrals  appearing  in  the  total  energy  equation  as 


6 


(3.77) 


z 


0 


6 


E 


(3.78) 


z 


0 
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The  above  equations  are  the  defining  equations  for  the  structure  functions 
Sj^d)  ,  s£b)  which  are  the  structure  functions  for  turbulent  decay  and  buoyant 
production/destruction  of  the  total  energy.  The  quantity  l  is  the  character¬ 
istic  turbulent  scale  size  with  the  boundary  layer.  This  quantity  is  governed 
by  the  integral  form  of  the  turbulent  scale  equation  which  is  described  in 
Section  3.7. 


In  the  usual  fasion  we  define  60  as  that  elevation  at  which  the 
viscous  turbulent  levels  have  fallen  to  an  arbitrarily  small  fraction  of  the 
maximum  values  within  the  boundary  layer.  We  choose  to  relate  this  thickness 
60  directly  to  the  total  energy  thickness  5^  .  Thus,  we  relate  6^  through 
a  parameter  r  as 


S0  =  r6E  (3.79) 

6  =  r6r  + 

E  o 


The  total  energy  integral  equation  now  becomes  an  equation  for  60  by 
eliminating  6^  through  Eq.  (3.79): 


9  (  r  r  \  _  (  d  ) 

9t  ^ErV  "  UE 


★ 

E  (6  - 

r  o 


«J  -  lo^'E 


4b>. 


r^o 


(3.80) 


In  Eq.  (3.80),  the  energy  decay  rate 


is  given  by 


O) 


(d)  _  ,(d) 
E  “  aE 


qr/(4&) 


(3.81) 


while  the  buoyant  production/decay  rate  is  given  by 


4b)  *  SEb>  I“*|3^kL  Er> 


(3.82) 


The  equilibrium  boundary  layer  thickness  6q  is  given  by 


*  4r*(4u„  +  v»VJ 

6o  - 


SEdVr 


(3.83) 
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We  see  that  in 
relaxes  to  the 


the  neutral  case  (L  °°) 
equilibrium  thickness  6* 


,  the 
on  a 


total  boundary  layer  thi^kjjiess 


time  scale  given  by 


thick 

4<i)- 
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3.3.7  Average  Scale  and  Integral  Scale  Equation 


It  will  be  noted  that  in  the  expression  of  the  various  integrals  of  the 
total  energy  equation  in  terms  of  non-dimensional  structure  functions,  it 
was  necessary  to  utilize  the  average  scale  £  .  We  now  take  up  the  structure 
function  transformation  of  the  integrated  scale  equation  (3.23)  which  governs 
£  .  We  define  the  structure  functions  for  the  turbulent  destruction  of 

seal e  as 


6 


z 


0 


(3.84) 


The  structure  function  for  turbulent  production  of  scale  is  defined 

by  the  statement 


6 


(3.85) 


z 


0 


The  structure  function  for  buoyant  production/destruction  of  scale  is  defined 
by 


6 


=9-  wu  dz  =  =  S 


o 


(3.86) 


z 
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The  integral  form  of  the  scale  equation  may  then  be  expressed  as 
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6  9t  <Wo> 
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(u*U  +  v*V  )  „  /n\ 

_  sr'  — —  4-+  s{p)q 


(d) 


qr  \\ 


l  T 


-  0.8  S 


(b)  ^ 


£  3 
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A 

L 


(3.87) 


We  may  express  the  above  result  as 


A  A  <tso>  £  -  "ib)t 


(3.88) 


where  -jS  the  scale  decay  rate  given  by 


(d>  _  .id)  (u*u»  +  v*vJ 

w£  "  b£  .  ifl  i 

\&0W 


(3.89) 


and  w^b)  is  the  buoyant  production/decay  rate  for  scale: 


w^b)  =  0.8  s[b)  — 


<q^.L 


The  equilibrium  scale  £  is  given  by 


(3.90) 


* 

£ 


(P)  2 


£ 

TdT 


q  \\ 


(uk 


v*V 

X  o 


(3.91) 


Nominal  values  for  these  structure  functions  are  s[  =  0.2,  = 

0.0125,  Sj£b)  =  0.5.  The  behavior  of  the  scale  in  the  neutral  and  stable 
limits  is  of  particular  interest.  In  the  neutral  limit  the  scale  will  take 
the  equilibrium  value  indicated  by  Eq.  (3.91).  For  the  above  nominal  values 
of  the  structure  functions  we  find  for  a  steady  state,  neutral  PBL  (with 
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(3.92) 


£*  =  0.06  6 
*  o 


in  good  agreement  with  the  characteristic  scale  predicted  in  Reference  24. 
In  the  strongly  stable  limit  (0  <  L  «  60)  the  scale  will  tend  to  an 
equilibrium  value  given  by  the  balance  between  the  last  two  terms  of 
eq.  (3.87) 


i 


s‘b)(0.8)|I, 


(3.93) 


which  yields  l  =  0.16  L  for  the  nominal  structure  function  values.  This 
result  is  in  agreement  with  the  appropriate  limit  on  £  of  ~  0.2  L  for 
strongly  stable  boundary  layers. 


3.3.8  Linearized  Forms  of  the  Integral  Equations 


We  may  express  the  integral  equations  for  the  mean  flow,  Eqs.  (3.20  - 
3.22)  in  terms  of  the  thickness  and  flux  vectors  6  ,  c  defined  in 

Section  3.4  as  a  a 


d6 


dt 


r  *  Itjjc  +  fZ  06q  -  B  060 
r  1  oo'  a  aB  B  aB  6 


(3.94) 


where  is  the  partial  anti-symmetric  matrix: 


Jap 


/  0  1  0 

-10  0 

\  0  0  0 


(3.95) 


and  is  the  matrix  of  forcing  functions: 


/f'"1  0  0  \ 


aB 


f(oo)  0 

I  y  U 


\°  0  fe”V 


(3.96) 


where 


f(°°)  =  f(°°)  _  _1 
TU  TV 


dt 


(3.97) 


/  \  t  d0 

F(°°)  _  J_  L 

6  ■  0r  dt 


(3.98) 


The  linear  expansion  of  ca  about  some  state  ca(0)  ,  Sa(0)  is 
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Ca  "  ca<°>  -  ^  WV0)) 


W^ere  ^a3  1S  t*ie  differential  transfer  matrix  defined  by  Eq.  (3.65).  Sub¬ 
stituting  this  form  into  Eq.  (3.94)  we  obtain 


d6  .  - 

-£■  =  |tJ  |c  (0)+  A  060  +  6_1  |UjC  r6  (0) 
dt  1  oo1  a  ag  g  o  ^  otp  p 


(3.99) 


where  Aag  is  the  fundamental  matrix  of  the  system: 


A  =  fZ  0  -  6_1 ifilC  a  -  B  , 
aB  a(3  o  1  001  a(3  aB 


(3.100) 


If  the  state  5a(0)  is  an  equilibrium  state,  then  the  equation  governing 
perturbations  6^  h  6a  -  Sa(0)  about  this  state  is 


d6’ 

—  =  A  6' 

dt  aB  p 


(3.101 ) 


We  thus  see  that  the  matrix  A^g  contains  all  the  fundamental  linear  response 
modes  of  the  system  including  tne  rotational  (Coriolis)  modes  and  the  Brunt 
Vaisaila  modes.  Two  of  the  eigenvalues  of  Aag  may  be  identified  with  the 
rotational  modes  while  the  third  eigenvalue  may  be  identified  with  the  Brunt 
Vaisaila  mode. 


The  integral  equations  for  the  boundary  layer  thickness  SQ  and  the 
scale  £  ,  Eqs.  (3.80)  and  (3.89)  are  similarly  in  a  form  in  which  linearization 
is  readily  applied.  We  point  out  the  linear  forms  of  these  equations  because 
they  form  the  basis  of  the  computational  solution  technique.  We  do  not  finite 
difference  the  equations  for  £„  ,  60  ,  £  .  Rather  we  utilize  the  linearized, 
constant  coefficient  forms  of  these  equations,  assuming  their  validity  over  a 
small  time  interval  At  connecting  the  two  states  at  t]  ,  t2  in  terms  of 
analytical  solutions  to  the  linear  forms  with  the  coefficients  held  fixed  at 
their  values  at  the  time  t-|  .  After  the  solution  is  obtained  at  time  t£  , 
the  coefficients  are  re-evaluated  and  the  process  is  repeated  for  the  next 
1 evel  tg  =  t2  +  At  ,  etc . 
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3.3.9  Coupling  of  the  Integral  and  Differential  Systems 

As  noted  previously,  the  integral  equations  for  the  variables 
V  <5y,  6q,  6,  £  and  the  fluxes  u* ,  v*,  0*  form  a  closed  dynamical 

system  if  the  structure  functions  S^,  S^\  S^,  Sr^,  S„  , 

si  '  are  known.  Because  of  the  manner  in  which  these  structure  functions 
are  defined,  they  are  only  sensitive  to  the  integrals  over  detailed 
profile  shape  and  are  not  sensitive  to  the  characteristic  magnitudes 
of  the  variables.  Thus,  one  can  specify  the  structure  functions  as 
pure  non-dimensional  numbers  and  obtain  reasonably  good  solutions 
for  the  integral  and  surface  parameters.  For  example,  for  a  neutral 
PBL  with  unsteady  forcing  which  is  not  too  rapid  compared  to  the  rotation 

rate  f,  the  values  =  0.8,  =  3.0,  =  0.15,  =  0.0125, 

=  0.2  r  =  0.35  will  yield  satisfactory  results  with  =  5I_C. 

For  flows  with  strong  forcing,  more  accurate  results  are  to  be 
obtained  by  introducing  a  finite  difference  system  in  the  "outer"  domain 
<5  <  z  <  6  wherein  the  general,  second  order  closure  system  of  equations 

is  solved. 

This  finite  difference  system  requires  boundary  data  at  the  surface 
z  =  6S  as  well  as  the  free  stream  conditions  at  z  =  6.  Finite  difference 
derivative  boundary  conditions  which  would  otherwise  be  required  at  the 
surface  layer  "edge"  z  =  are  dispensed  with;  instead,  the  boundary 
data  for  u(6z,  t),  0 (6S  >  t),  ...  uw(6s,  t),  ...  are  taken  from  the  integrally 
determined  surface  layer  conditions  and  flux  laws.  Hence,  the  surface 
boundary  condition  for  the  differential  equation  set  are  fixed  by  the 
integral  parameters. 

The  coupling  back  upon  the  integral  constraint  equations  from  the 
finite  difference  generated  profiles  is  through  the  non-dimensional 
structure  functions  which  involve  integrals  over  the  profiles  of  the 
velocity,  temperature,  and  turbulence  fields.  In  the  limit  in  which 
no  finite  difference  points  are  used,  i.e.,  the  method  is  purely  integral, 
these  structure  functions  may  be  specified  as  pure  numbers.  Hence,  at 
appropriate  time  levels  in  the  course  of  evolution  of  the  integral 
equations,  the  structure  functions  are  up-dated  by  explicit  calculation 
of  the  integrals  over  the  profiles  generated  by  the  evolution  of  the 
differential  equation  set. 

The  finite  difference  system  is  set  up  on  a  dynamically  moving  grid 
whose  first  point  is  located  at  z  =  6s(t)  and  whose  uppermost  (top)  point 
is  located  at  z  =  6(t).  The  motion  of  the  grid  is  thus  fixed  by  6(t) 
which  is  in  turn  determined  by  the  integral  equation  of  the  total  kinetic 
energy.  The  basic  flow  of  information  is  shown  schematically  in  Figure  3.1. 
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Figure  3.1  -  Information  flow  in  hybrid  integral -differentiation 
description  of  the  PBL. 
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3.4  Illustrations 


3.4.1  PBL  Subjected  to  Unsteady  Forcing 

As  a  basis  for  illustration  we  consider  a  particular  unsteady  problem  for 
the  neutral  PBL.  These  problems  will  involve  ramp  transitions  of  the  inviscid 
region  velocity  from  one  steady  state  level  to  another.  If  U^O)  is  the 
initial  steady  state  and  U^U)  is  the  final  steady  state,  the  ramp  transition 
is  defined  by 


UJO) 


t  <  tf 


U„(t)  -  <U„(0)  +  M1)  '  M°--  f(t-t„)  to  <  t  <  t0+  tf-1  (3.102) 


u.U) 


Tf"1  <  t 


The  transient  is  specified  by  the  parameters  U^O),  U^U),  x.  The  general 
problem  is:  A  boundary  layer  in  steady  state  corresponding  to  an  inviscid 
velocity  U^O)  is  subjected  to  a  linearly  increasing  inviscid  velocity  over  a 
time  period  if  until  it  reaches  a  value  U00(l)  at  which  point  the  velocity 
remains  fixed  at  the  new  steady  state  value  U^l).  Determine  the  motion 
within  the  PBL.  As  such,  this  problem  allows  us  to  study  the  neutral  PBL 
subject  to  both  Coriolis  effects  and  unsteady  forcing. 

The  geostrophic  conditons  (Ug,  Vg)  are  established  so  that  for  all  time. 


U 


oo 


ujt) 


V  =  0 


oo 


(3.103) 


We  select  conditions  exhibited  in  Table  3.4.1.  It  should  be  observed 
that  these  conditions  describe  a  very  severe  transient  in  that  the  forcing 
frequency  f^u)  [Eq.  (3.97)]  is  initially  ten  times  greater  than  the  Coriolis 
frequency  f. 
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Table  3.4.1 


Conditions  for  Neutral  PBL  Subjected  to  Ramp  Transition 
U«,(0)  =  1  m/sec 
U0O(l)  =  10  m/sec 

T  =  1 

f  =  10-1*  sec-1 


The  solutions  for  this  problem  are  determined  in  three  different  ways.  In  the 
first,  we  utilize  the  standard  A.R.A.P.  second-order  closure  theory 
implemented  in  a  fully  finite-difference  procedure  utilizing  of  the  order  of 
40  grid  levels.  The  PBL  response  (as  reflected  in  surface  RMS  turbulence  and 
cross-isobaric  angle)  computed  in  this  manner  is  shown  in  Figure  3.2.  In  the 
coordinates  utilized,  the  Coriolis  period  is  equal  to  2u.  It  can  be  seen  that 
q0  makes  an  initially  rapid  transition  (including  an  overshoot  and  undershoot) 
during  the  period  of  acceleration  of  the  outer  flow  and  then  oscillates  about 
its  new  level  with  the  oscillation  slowly  damping.  The  angle  particularly 
evidences  higher  harmonics  of  f.  These  result  from  the  nonequilibrium 
rotational  wave  modes  which  are  present  in  the  turbulence  equations. 

In  Figure  3.3  we  show  the  response  to  the  same  problem  as  computed 
in  the  second  manner:  the  hybrid  procedure  with  no  finite  difference 
levels,  i.e.,  the  procedure  in  pure  integral  form.  The  integral  model 
accurately  exhibits  the  initial  overshoot  in  qo  (but  fails  to  give  an 
undershoot),  and  then  yields  a  similar  decline  over  the  period  of 
acceleration.  Because  it  does  not  contain  dynamic  equations  for  the 
full  Reynolds  stresses,  the  pure  integral  form  of  the  hybrid  model  does 
not  yield  the  higher  harmonics  of  f,  but  only  oscillates  at  the  fundamental. 
The  average  angle  response  follows  the  full  finite  difference  solution 
well  but  climbs  more  rapidly  to  the  peak  value  following  the  period  of 
outer  flow  acceleration. 

In  Figure  3.4  we  show  the  response  to  the  same  problem  as  computed 
in  the  third  manner:  the  hybrid  procedure  with  5  finite  difference  levels. 

The  response  is  equally  well  predicted  with  a  tendency  for  the  first 
harmonic  (but  not  a  second)  to  appear  in  the  angle  response.  We  remark 
that  the  hybrid  procedure  with  only  5  grid  levels  is  somewhat  sensitive 
to  numerical  instabilities  in  the  following  sense.  If  the  profile 
computation  over  5  points  develops  any  significant  errors,  the  degrading 
influence  on  the  structure  functions  can  feed  back  through  the  surface 
conditions  and  further  degrade  the  profile  structure  near  the  surface. 

Nonetheless,  the  hybrid  procedure  executes  at  least  8  times  faster 
than  the  full  finite  difference  procedure  and  as  much  as  40  times  faster 
for  the  purely  integral  version.  The  good  quality  of  these  results, 
given  the  vast  decrease  in  computer  resources  spent,  seems  highly 
worthwni le. 
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Standard  A.R.A.P.  Finite  Difference  Model 


I  75 


Finite  Difference  Levels:  30-40 
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Figure  3.2(a).  Surface  RMS  turbulence  response  to  unsteady  forcing  of  the  PBL. 

Period  of  outerflow  acceleration  is  from  f t  =  0  to  f t  =  1 . 
Standard  ARAP  finite  difference  model. 


Angle  (  Deg. ) 


Hybrid  Model 

Finite  Difference  Levels-  0 


Time  ( tf ) 

Figure  3.3(b).  Cross  isobaric  angle  for  conditions  of  Figure  3.3(a). 
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Figure  3.4(b).  Cross  isobaric  angle  for  conditions  of  Figure  3.4(a). 


3.4.2  Stable  PBL  Illustrations 


As  a  second  illustration,  we  consider  the  computation  of  stable  PBL 
response.  Our  illustration  for  this  case  will  be  the  archetypal  problem 
of  a  constant  surface  cooling  rate  in  which  a  quasi -steady  boundary  layer 
is  established  (References  25-27).  We  begin  with  an  equilibrium  neutral 
PBL  at  time  t  =  0  at  which  time  a  constant  surface  cooling  rate  is 
applied.  After  a  period  of  several  Coriolis  periods,  the  quasi -steady 
state  is  established.  The  computation  illustrated  here  is  for  the 
purely  integral  version  of  the  hybrid  procedure.  The  conditions 
for  the  illustration  are  presented  in  Table  3.4.2. 


Table  3.4.2 

Conditions  for  Constant  Surface  Cooling 
Rate  Stable  PBL  Illustration 
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The  evolution  of  the  boundary  layer  for  these  conditions  is 
shown  in  Figure  3.5.  After  a  transient  of  approximately  2  Coriolis 
periods  the  characteristic  boundary  layer  parameters  approach  quasi  - 
steady  values.  These  quasi-steady  values  are  consistent  with  those 
predicted  in  Reference  25.  The  hybrid  model  indicates  an  initial 
undershoot  in  cross-i sobaric  angle  of  -15°  before  evolving  to  a  steady 
state  values  of  approximately  50°.  This  steady  state  value  is 
about  10°  larger  than  that  predicted  in  Reference  25;  no  attempt  has 
been  made  to  fine-tune  the  values  of  the  structure  functions  in  these 
illustrations  to  effect  more  exact  comparisons.  The  value  of  the 

Zilitinkevich  parameter  d  =<VV/|u*lL  is  shown  in  Figure  3.5(e)  and 

its  value  of  0.44  at  quasi-steady  state  is  closely  consistent  with 
that  of  Reference  24  and  25. 
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surface  cooling  rate  corresponding  to  the  condition  of  Table  3.4.2. 


I  u*l  /  I  Ucol 


Angle  (Deg.) 


Figure  3.5(c).  Cross-isobaric  angle  response  corresponding  to  Figure  3.5(a). 
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4.  THE  INTERACTION  OF  TURBULENCE  WITH  PRECIPITATION: 
FORMULATION  OF  A  PRECIPITATION  MODEL 
FOR  THE  PLANETARY  BOUNDARY  LAYER 


4.1  Introduction 

Many  of  the  atmospheric  Planetary  boundary  layer  (PBL)  processes  and 
episodes  which  are  successfully  discrirable  in  terms  of  a  second  order 
closure  theory  of  turbulence  such  as  that  in  use  at  ARAP  (Reference  10) 
involve  the  transport,  condensation,  and  evaporation  of  water.  It  is 
important  that  such  models  correctly  describe  the  overall  balance  of 
water  by  the  various  evaporation,  condensation,  diffusion,  and  advection 
mechanisims.  In  some  cases  of  interest,  precipitation  of  cloud  water  in 
the  form  of  rain  or  drizzle  to  the  surface  is  an  important  process  controlling 
the  balance  of  water.  In  addition,  precipitation  is  of  interest  in  its 
own  right;  for  it  is  desirable  to  predict  the  likelihood  and  magnitude  of 
precipitation  drizzle  or  rain  events.  In  the  present  work  we  therefore 
present  a  discussion  and  model  of  precipitation  for  use  in  atmospheric  PBL 
models.  Because  we  restrict  attention  to  the  PBL,  we  consider  only  warm 
(non-freezing)  precipitation  and  cloud  droplet  growth  processes. 

Following  the  stage  of  condensation  growth  in  the  early  period  of 
cloud  formation  (droplet  radii  r  <  10  ym),  further  growth  of  the  cloud 
droplets  to  reach  precipitation  size  is  generally  believed  to  occur  by 
colli  si onal  coalescence  of  droplets.  Gravitational  sedimentation  has 
received  virtually  the  exclusive  attention  of  theorists  as  the  colli  si onal 
coalescence  mechanism  of  atmospheric  clouds  (References  28-30).  There 
seems  little  doubt  that  the  collisional  coalescence  of  drops  of  different 
size  is  an  important  droplet  growth  mechanism  at  some  stages  of  cloud 
evolution.  On  the  other  hand,  in  the  early  stages  of  growth  (1  <  r  <  50  ym) 
this  mechanism  possesses  certain  limitations.  Two  of  these  limitations 
are  (1)  the  inherent  requirement  of  differential  size  for  a  non-zero 
collision  rate,  and  (2)  the  sharply  diminished  collision  efficiencies 
which  result  for  the  vanishing  relative  Reynolds  numbers  of  differential 
sedimentation  when  both  collision  partners  approach  the  same  size.  These 
two  limitations  when  viewed  in  the  light  of  a  further  result  from  classical 
condensation  theory  -  namely  the  narrowing  of  the  droplet  spectrum  into  a 
single  size  during  condensation  growth  -  suggest  that  collisional  growth 
mechanisms  other  than  gravitational  sedimentation  may  play  an  important 
role  in  the  initial  growth  stage  of  clouds  into  precipitation  size  drops. 

In  particular,  atmospheric  turbulence  may  play  an  important  and  direct  role 
in  the  evolution  of  the  cloud  drop  spectrum. 

We  propose,  therefore,  to  include  in  the  collision  mechanism  of  our 
precipitation  model  the  effects  of  atmospheric  turbulence  in  addition  to 
classical  gravitational  sedimentation.  There  appear  to  be  at  least  two 
unique  ways  in  which  turbulence  affects  the  evolution  of  the  cloud  drop 
spectrum.  The  first  is  in  the  dynamics  of  condensation.  Although  classical 
condensation  theory  predicts  a  narrowing  spectrum,  (forcing  a  single  size 
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cloud  droplet)  it  appears  that  when  condensation  is  considered  in  a  turbulent 
environment,  the  spectrum  may  be  broadened  (References  31,  32).  The  second 
role  is  in  the  collisional  growth  stage  in  which  the  turbulent  velocity 
field  provides  not  one  but  several  mechanisms  for  the  collisions  of  drops. 

In  Section  4.2  we  review  the  stages  of  condensation  and  growth  processes 
for  clouds  pointing  out  the  regimes  in  which  turbulence  can  be  of  importance. 

In  Section  4.3  we  discuss  in  detail  the  turbulence  collision  mechanisms  and 
develop  the  collision  kernel  for  the  processes  of  both  turbulence  induced 
collisions  and  gravitational  sedimentation.  In  section  4  we  present  a  two- 
group  precipitation  model  based  upon  the  growth  mechanisms  discussed  in  Section 
4.3  which  we  term  the  Cloud-Precipitation  (CP)  model.  In  Section  4.5_ 
we  present  some  illustration  calculations  of  the  CP  model  for  the  infinite, 

homogenous  cloud. 

4.2  Condensation,  Evolution  of  the  Cloud  Droplet  Spectrum,  and  Precipitation 

To  set  the  stage  for  the  model  we  propose,  we  first  review  the  various 
processes  which  take  place  from  the  onset  of  a  water  mixing  ratio  in  excess  of 
the  saturation  value  to  the  final  stage  (if  it  occurs  in  the  time  scale  of  a 
particular  problem)  in  which  drops  precipitate  to  the  surface.  We  point  out 
the  role  of  turbulence  in  certain  of  these  stages  of  development.  The  stages 
of  drop  evolution  may  be  defined  in  the  following  scheme: 


(1)  The  Nuclei  Activation  Stage  _2 

drop  (particle)  size  =  10  pm  <  r  <  1  pm 

time  scale  =  1  sec.  or  less 


(2)  The  Condensation  Growth  Stage 

drop  size 
time  scale 

(3)  The  Collisional  Growth  Stage 

drop  size 
time  scale 


1  pm  <  r  <  10  pm 
1-100  sec. 


10  pm  <  r  <  10 3  pm 
>  100  sec.  but  highly  variable 


(4)  The  Sedimentation  Stage  w 

drop  size  =  100  pm  <  r  <  103  pm 

time  scale  =  1 o ~ 3 -1 0  hours 

Once  the  water  mixing  ratio  exceeds  the  local  saturation  value,  nuclei 
must  be  activated  before  water  may  condense  in  realistic  time  scales.  Follow¬ 
ing  nuclei  activation,  drops  grow  by  the  direct  condensation  of  vapor  and  in¬ 
teract  negligibly  via  collisions.  In  most  atmospheric  situations,  the  liquid 
water  formed  by  the  overall  amount  of  excess  saturation  and  the  number  of  nu¬ 
clei  available  and  activated  results  in  a  cloud  with  drop  number  densities 
ranging  from  108-109  nf3  and  radii  ranging  from  1  to  10  pm.  It  should  be  noted 
that  classical  condensation  theory  predicts  a  spectrum  evolution  through  the 
stages  of  nuclei  activation  and  condensation  growth  which  is  progressively 
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narrowing  (Reference  33)  ultimately  yielding  a  cloud  in  which  all  particles 
are  the  same  size.  Thus,  in  classical  condensation  theory,  size  differential 
in  the  cloud  spectrum  can  only  result  from  the  residue  of  size  differential 
of  nuclei.  A  limited  amount  of  work  has  been  devoted  to  the  examination 
of  droplet  evolution  in  the  nuclei  activation  and  condensation  growth  stages 
in  a  turbulent  environment  (References  31,32).  Although  the  "turbulence" 
models  used  in  these  investigations  are  seriously  over  simplified,  it  does 
appear  that  turbulent  condensation  theory  provides  significant  broadening 
of  the  droplet  spectrum  which  can  dominate  over  the  "natural  narrowing" 
of  classical  condensation  theory.  We  have  made  some  preliminary  investigations 
of  more  rigorous  turbulent,  second  order  closure  versions  of  the  droplet 
kinetic  equations  and  find  both  broadening  and  narrowing  effects  which  can 
result  from  the  turbulent  correlations  w"' w1  ,W‘ e 1  where  w  is  the  vertical 

velocity,  e  the  potential  temperature,  and  an  overbar  denotes  a  turbulent 
ensemble  average.  We  shall  report  on  these  investigations  in  subsequent 
publicati ons. 

For  cloud  droplets  to  grow  significantly  beyond  the  range  of  radii 
r  <  10  ym  up  to  precipitation  sizes,  collisional  coalescence  processes 
must  be  operative.  The  time  scale  for  these  col  1 i sional  processes  are  highly 
vari able  when  turbulence  induced  collisions  are  included  in  addition  to  the 
gravitational  sedimentation  collisions  of  classical  coalescence  theory.  We 
find  this  result  (which  will  be  discussed  in  Section  3)  consistent  with  the 
highly  variable  nature  of  natural  clouds  to  grow  to  rain  or  drizzle  size  drops. 
A  critical  feature  of  the  collisional  growth  process  is  the  creation  of  a 
small  number  of  drops  much  larger  than  the  average.  This  long  tail  effect 
in  the  distribution  function  is  the  result  of  the  increasing  collision 
cross-section  of  large  drops. 

The  final  stages  of  drop  evolution  occur  when  drops  have  grown  signi¬ 
ficantly  large  enough  to  develop  a  significant  precipitation  velocity.  These 
preci pi  table  drops  then  leave  the  cloud  and  progress  to  the  surface  where 
they  then  leave  the  atmosphere. 

The  precipitation  model  we  shall  describe  in  section  4  in  the  present 
work  treats  the  first  two  stages  of  cloud  evolution  described  above  in  very 
simple  parametric  fashion.  The  collisional  growth  and  sedimentation  processes, 
however,  will  be  treated  in  some  mechanistic  detail  that  includes  the  effects 
of  turbulence  upon  collisional  growth. 

4.3  Turbulence  and  the  Coll isional/Coal escence  Process 

We  first  consider  the  conceptual  picture  of  atmospheric  turbulence  and 
its  influence  on  the  collisions  of  liquid  drops  embedded  in  such  an  environ¬ 
ment.  Atmospheric  turbulence  consists  of  the  random  motion  of  eddy  struc¬ 
tures  ranging  from  the  largest  energy  containing  scales  to  the  dissipation 
scales  where  molecular  viscosity  comes  into  play.  The  largest  scale  is  of  the 
order  of  the  largest  characteristic  macro-length  L  (e.g.  the  PBL  depth,  ter¬ 
rain  dimension,  etc.  )  while  the  smallest  scale  is  of  the  order  of  the  Taylor 
micro-scale  A.0  defined  in  terms  of  the  turbulence  energy  dissipation  e  =  o3/A 
as  Xo  =  (n  /e)%  where  p  is  the  molecular  kinematic  viscosity.  Similarly^  a 
micro-time  characterizing  the  time  scale  of  fluctuations  of  the  dissipation  can 
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be  defined  as  t0  =(n/e)%  and  a  micro-acceleration  cnaracteri zing  tne  accel¬ 
eration  of  the  flow  field  within  the  eddie  as  a0  =  a0/t02.  The  magnitudes  of 
these  quantities  for  the  range  of  turbulence  dissipation  rates  encountered  in 
the  atmosphere  are  present  in  table  4.1. 


Table  4.1 

Turbulence  Length,  Time,  and  Acceleration  Scales  in  the  Atmosphere 

for  Dissipation  Scale  Eddies 


£ 

Ao 

TO 

ao/q 

(m2/sec. 3) 

(  ym  ) 

(  ms  ) 

_ 

0.001 

1510 

130 

0.Q09 

0.01 

846 

42 

0.05 

0.10 

476 

13 

0.28 

1.0 

268 

4.2 

1.58 

10.0 

151 

1.3 

8.9 

Since  drops  of  radius  r  will  generally  satisfy  the  condition  r  «  L, 
the  nature  of  turbulence  induced  collisions  will  turn  first  on  the  question  of 
whether  r  >  A0  or  r  <  Ao . 

It  can  be  seen  that  even  for  the  highest  dissipation  rates,  only  drops 
greater  than  about  150  ynr  would  be  larger  than  the  dissipation  scale  eddies. 
We  conclude,  therefore,  that  under  most  circumstances  (and  parti cul ary  for 
drop  sizes  in  the  crital  range  1  <  r  <  50  ym  )  a  cloud  droplet  will  execute 
its  collisional  dynamics  within  a  dissipation  scale  eddy.  The  precise  nature 
of  the  flow  field  within  a  dissipation  scale  eddy  is  not  clearly  understood 
at  present;  however,  such  flow  fields  must,  of  necessity,  be  characterized 
by  a  high  shear  rate. 

The  average  shear  rates  in  the  dissipation  scale  ed.dy  can  be  related 
to  the  turbulence  dissipation  as  (Reference  34) 

S  s  (e/n) "  =  Vto 

The  first  effect  of  turbulence  upon  the  collision  rate  of  cloud  drops 
is  thus  to  place  them  in  a  shear  flow  of  average  shearing  rate  S.  Thus,  two 
drops  of  radii  rx ,  rx  lying  within  a  collision  cylinder  will  possess  a  re¬ 
lative  velocity  with  respect  to  one  another  of  magnitude 

A =  (r,  +  r2)  s 

This  collisional  relative  velocity  mechanism,  in  contrast  to  that  of  gravi¬ 
tational  sedimentation,  does  not  require  a  size  difference  between  the  col¬ 
lisional  partners. 

The  second  effect  of  turbulence  is  to  create  an  acceleration  field  a0 
for  the  flow  field  of  the  drops  in  addition  to  that  of  gravity.  Hence,  dif¬ 
ferential  size  relative  motion  will  be  enhanced  by  the  presence  of  the  tur¬ 
bulent  field. 
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The  third  effect  of  turbulence  upon  colli  si onal  encounter  is  for  the 
drops  with  r  >  \q.  This  regime  is  the  most  complex  involving  the  com¬ 
plication  of  the  flow  field  of  eddy  scales  larger  than  the  Taylor  scale. 

A  simulation  of  the  collision  dynamics  which  may  in  some  respects  model 
this  regime  with  an  arbitrarily  varying  background  velocity  field  but 
with  turbulence  induced  shearing  neglected  Almeida  (Reference  35)  indicates, 
as  one  would  expect,  the  enhancement  of  colli  si onal  efficiencies  of 
such  drops  due  to  the  agitation  of  the  background  field.  However, 
the  collision  efficiencies  of  such  large  drops  are  already  much  larger 
than  the  minimal  levels  characteristic  of  drops  in  the  1-50  ym  range 
even  in  the  absence  of  turbulence.  In  addition,  the  critical  range 
for  growth  of  1-50  ym  will  be  such  that  drops  turbulently  collide 
primarily  by  the  shearing  mechanism.  We  thus  disregard  the  turbulence 
effects  on  drops  with  r  >  Ao  and  include  the  effects  of  turbulence 
on  the  shearing  rate  and  acceleration  of  the  flow  field  surrounding 
drops  within  dissipation  scale  eddies. 

Let  us  now  consider  the  formulation  of  the  collision  kernel  for  the 
col  1 i sional  encounters  of  such  drops.  The  total  collision  rate  per  unit 
volume  between  two  populations  of  drops  of  radius  n  and  number  density 
n i  and  radius  r2  and  number  density  n2  may  be  expressed  as 

N  =  n!  n2  V0  (rlsr2) 

where  v0  (rx,  r2)  is  the  collision  kernel.  Five  basic  processes  con¬ 
tribute  to  the  collision  kernel  in  the  atmosphere.  These  are 

1.  Turbulent  Shearing 

2.  Turbulent  Accelerations 

3.  Gravitational  Sedimentation 

4.  Brownian  Motion 

5.  Electrostatic  Attraction 

In  what  follows,  we  restrict  attention  to  the  first  three  processes.  Brownian 
motion  is  only  expected  to  be  important  for  particle  sizes  much  smaller  than 
the  average  cloud  drop.  No  attempt  is  made  to  estimate  the  influence  of  elec¬ 
trostatic  attraction. 

The  collision  kernel  v0  may  be  expressed  as  (Reference  36) 


Vfi2  =  V!  2  +  V22  +  V32 


Where  'o1 ,  v2  s  v3  and  the  contributions  of  turbulence  induced  shearing, 
turbulence  induced  accelerations,  and  gravitational  sedimentation  resDPc- 
tively  .  These  are  expressed  as 


v,  =  AV12(s)  Ej  (rx,  r2h  rf2 

v2  =  /TT3^AV12(s)  F 2  (rl5  r2)  tt  rzlz 

v3  =  AV12(8)  E3  (rx,  r2 )ir  r\2 


69 


In  the  above  AVi^  is  the  effective  relative  velocity  of  the  drops  due  to 
the  shearing  motion  of  the  fluid  expressed  in  terms  of  the  shearing  rate  S: 

A vx|s)  =  r!2  s 

where  S  is  given  in  terms  of  the  turbulent  dissipation  rate.  The  collision 
cylinder  radius  is  r12  =  (n  +  r2).  The  relative  velocity  A \lAz)  is  the 
difference  between  the  terminal  gravitational  sedimentation  velocities  of 
the  two  particles: 

iVl2(g)=  IV,  -v2| 

The  form  of  V.^  (r-j)  depends  upon  the  flow  regime  of  the  particles.  In  the 
Stokes  range 

Vi  (ri>  =(ft  I  ?ri2 

where  pQ  are  the  densities  of  drop  liquid  and  surrounding  fluid  respec¬ 
tively  and  n  is  the  kinematic  viscosity  of  the  surrounding  fluid.  The  quan¬ 
tities  E1,  E2,  E3  are  the  collision  efficiencies  for  each  of  the  processes 

respecti vely.  These  are  defined  in  terms  of^the  cross-sectional  area  Q  per¬ 
pendicular  to  the  relative  velocity  vector  AV12  within  which  the  centers  of 
the  drops  must  lie  if  they  are  to  collide  compared  to  the  geometrical  hard  sphere 
collision  cross-section  of  the  two  drops: 


Although  it  is  not  indicated  functionally,  the  collision  efficiencies 
Ei  are  functions  of  the  radius  ratio  of  the  colliding  drops  and  the  relative 
Reynolds  number.  For  gravitational  sedimentation,  this  number  depends  pure¬ 
ly  upon  the  radii  rls  r2.  For  shearing  motion,  however,  the  efficiency 
depends  upon  the  shearing  rate  as  well  as  the  radii  rj  r2-  The  collision 
efficiency  E3  in  the  absence  of  turbulence  is  summarized  in  Reference  37. 
There  are,  as  yet,  no  reliable  calculations  or  measurements  of  the  efficiencies 
E19  E2  (or  correspondingly  an  overall  efficiency  which  depends  upon  shearing 
rate  )  although  the  work  in  Reference  38  is  noteworthy. 

We  thus  may  summarize  the  three  colli sional  processes  of  interest  here 
by  noting  that  of  the  three,  shearing  collisions  are  the  only  ones  which 
are  operative  among  drops  of  equal  size;  hence  turbulence  provides  a  mecha¬ 
nism  (outside  of  the  Brownian  range)  of  coagulating  drops  of  equal  size  which 
is  otherwise  not  available  in  the  more  conventional  gravitational  sedimenta¬ 
tion  picture  of  collisional  coagulation. 
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4.4  The  Cloud  Precipitation  (CP)  Model 

4.4.1  The  Division  of  Liquid  Water  Into  Cloud  and  Precipitation  Groups 

As  described  in  section  2,  the  collisional  evolution  of  the  droplet 
spectrum  of  the  cloud  is  such  that  a  small  number  of  larger  droplets  are 
created  and  characterized  as  the  tail  of  the  cloud  distribution.  Because 
the  droplet  sedimentation  velocity  is  a  strong  function  of  particle  size,  it 
is  not  useful  to  characterize  the  precipitation  flux  as  an  average  over 
the  entire  liquid  water  distribution.  This  is  because  the  bulk  of  the  drop¬ 
lets  have  negligible  sedimentation  velocities.  Rather,  it  is  useful  to  de¬ 
fine  a  precipitation  group  as  those  particles  with  sedimentation  velocities 
greater  than  a  certain  minimum  value.  This  minimum  value  cannot  be  given 
by  the  cloud  micro  physical  process,  but  is  determined  by  the  overall 
macro-dynamics  of  the  problem  at  hand.  This  sedimentation  velocity  is 
selected  so  that  a  drop  will  fall  over  a  characteristic  macro-length  in  some 
characteristic  macro-time.  We  thus  select  (as  a  model  parameter  for  the 
precipitation  process)  a  sedimentation  velocity  (or  corresponding  parti¬ 
cle  size)  which  separates  cloud  droplets  (whose  contribution  to  the  sedi¬ 
mentation  flux  we  neglect)  from  drizzle  or  rain  drops  (which  constitute 
the  full  precipitation  flux.)  Let  us  designate  this  velocity  as  V*  and 
the  corresponding  particle  radius  as  r*.  Since  most  dynamical  events 
within  the  PBL  take  place  on  a  length  scale  of  the  order  of  10  m  or  less 
and  on  a  time  scale  of  the  order  of  1  hr-  we  find  the  minimum  precipita¬ 
tion  velocity  should  be  greater  than  or  equal  to  about  1  km/ hr  which,  cor¬ 
responds  to  the  sedimentation  velocity  of  a  particle  of  r*~  50  pm  in  still 

air.  The  total  liquid  water  is  thus  divided  into  two  groups:  a  cloud 
group  consisting  of  all  droplets  with  sizes  r  <  r*  and  a  precipitation 
group  with  sizes  r  >  r*-  Let  us  now  specify  the  various  collisional  and 
and  condensation  processes  which  take  place  between  these  two  groups. 

We  choose  not  to  describe  the  details  of  nuclei  activation  and  the 
dynamics  of  cloud  spectrum  formation.  These  processes  may  be  summarized 
in  terms  of  two  model  parameters:  the  average  cloud  droplet  radius  Rc 
and  the  non  dimensional  dispersion  of  the  cloud  spectrum  ac.  For  the 
present  model,  we  choose  to  consider  the  limit  ac  =  0  and  to  specify  Rc 
as  the  single  cloud-type  parameter  which  for  virtually  the  full  range  of 
cloud  types  lies  in  the  range  5  pm<:Rc<20  pm.  Thus,  given  the  total 
liquid  water  present  as  cloud,  the  cloud  droplet  number  density  nc  is 
implied  in  terms  of  the  cloud  droplet  average  radius  R^.  It  also  becomes 
clear  that  in  addition  to  the  usual  conservation  equations  for  total  liquid 
water,  two  additional  conservation  equations  are  required  to  determine  the 
precipitation  water  content  and  the  number  density  (or  average  size)  of 
precipitation  drops. 

With  the  cloud  droplet  variables  "c,  Rc,  so  determined,  three  col-, 
lisional  interaction  processes  and  an  evaporation  process,  then  emerge  which 
define  the  precipitation  drop  group  characteristics.  These  three  collisional 
processes  are  the  cloud-cloud  collisions,  cloud-precipitation  collisions, 
and  precipitation-precipitation  collisions.  Cloud-cloud  collisions  whose 
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coalescences  result  in  drops  with  radii  greater  than  r*  constitute  the 
cloud -  to -  cloud  precipitation  conversion  process.  Cloud-precipitation- 
precipitation  collisions  constitute  the  precipitation  aggregation  process. 
This  latter  process  is  only  important  in  situations  when  the  precipitation 
drop  number  density  is  very  large.  When  precipitation  drops  exist  in  other¬ 
wise  unsatuated  air,  evaporation  of  the  precipitation  drops  takes  place,  and 
we  term  this  process  the  precipitation  evaporation  process. 


4.4.2  Cloud  and  Precipitation  Variables  and  Conservation  Equations 

The  usual  mixing  ratios  are  given  in  terms  of  the  average  droplet  sizes 
Rc,  Rp  and  number  densities  nc,  np  for  cloud  and  precipitation  groups  res¬ 
pectively  as 


Hc  =  4/3  7iRc3  Pp/Poo  nc 
HP  =  4/3  Try  Po/P^  "p 
H*  =  Hc  +  Hp 
H  =  Hv  +  H£ 


(4.1) 


In  Eqs.  (4.1)  Pv  is  the  mass  density  of  water  in  the  vapor  phase, 
while  p*  is  the  mass  density  of  the  liquid  water,  and  p*.  is  the  mass 
density  of  the  mixture.  The  mixing  ratio  for  vapor  is  Hv>  that  of 
cloud  water  Hc,  and  that  of  precipitation  water  Hp.  The  total  liquid 
mixing  ratio  is  Hp  and  the  total  water  mixing  ratio  is  H. 

The  relationship  between  the  constituents  is  as  follows.  If  the  mix¬ 
ture  is  unsaturated  (H  <  Hs),  then 


H£  =  Hp 
Hy  =  H  -  H  '£ 


(4.2) 


In  these  statements  it  is  assumed  that  the  cloud  droplets  are  in  equili¬ 
brium  with  the  vapor;  the  precipitation  drops  need  not  be  in  equilibrium. 
If  the  mixture  is  saturated  (H  >  Hs)  then 


H  v  =  H  s 


Hc  =  H  -  Hs  -  Hp 
H£  =  Hc  +  Hp 


(4.3) 


For  the  general  case,  we  may  thus  write 
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(4.4) 


Hc  =  (H  -  Hs  -  Hp)^(H  -  H§) 

Hz  =  Hc  +  Up 

Hv  =  H  -  H£ 

Where  ^  (x)  is  the  Heaviside  function. 

It  will  be  noted  that  in  addition  to  the  saturation  mixing  ratio  Hs, 
the  water-air  mixture  system  possesses  two  degrees  of  freedom  as  we  have 
constructed  it  consisting  of  cloud  drops  and  rain  drops. _  If  H  and  HU 
are  specified,  all  the  water  species  variables  are  determined.  In  general 
we  shall  utilize  H  and  Hp  as  the  two  independent  variables  which  de¬ 
termine  the  various  water  species  variables.  It  can  be  readily  seen  that 
these  variables  are  sufficient  to  fix  the  values  of  Hv,  Hc  >  and 


The  precipitation  drop  number  density  conservation  statement  is 

Dn  „ 

_E.  +  3 —  (n  V  •) 

H4-  CW  V  p  pi  ' 


8xi 


Ncp  -  Npa 


(4.5) 


In  the  above  Vm  is  the  average  sedimentation  velocity  of  the  precipi¬ 
tation  drops,  i he  production  of  precipitation  drops  from  cloud  droplets 
by  the  cloud  to  precipitation  conversion  process  is  Ncp-  The  loss  of 
precipitation  drops  by  self-collisions  among  the  precipitation  drops  is 
Npa. 


The  conservation  equation  of  precipitation  water  is 


(hpV  ■ 


ft  +  H 
cp 


cc 


pe 


(4.6) 


The  production  rate  of . precipi tation  water  by  the  cloud  to  precipitation 
conversion  process  is  Hep*  The  production  of  precipitation  water  by  the 
cloud  collection  process  is  Hcc*  The  loss  of  precipitation  water  by 
evaporation  of  precipitation  drops  in  unsaturated  air  is  Hpe*  The  various 
rates  Ncp,  Npa,  Hep,  Hcc,  Hpe  are  described  in  subsequent  sections.  Equa¬ 
tions  (4.5)  and  (4.6)  provide  the  additional  dynamical  equations  which 
fix  the  properties  of  the  precipitation  group,  since  these  equations  deter¬ 
mine  np  and  Hp,  the  average  precipitation  drop  radius  Rp  is  determined 
the  third  of  Eqs.  (4.1).  The  conservation  equation  of  total  water  mixing 
ratio  is  modified  by  the  presence  of  a  precipitation  flux  and  becomes 


K  +  W  =  0  (4-7) 

4.4.3  The  Cloud  Conversion  to  Precipitation  (CP)  Process 

We  now  describe  the  CP  process  and  develop  expressions  for  the  CP  rates 
Ncp»  Hep*  Let  v  be  the  volume  of  any  given  drop  and  let  Vj  be  the  volume 
of  the  smallest  cloud  drop  under  consideration.  Then  m  =  v/vt  =  (yri) 
is  a  size  specification  parameter.  The  number  density  of  drops  of  size 
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m  we  denote  as 
tation  drops  np 


tim .  The  total  number  of  cloud  droplets 
are  then  given  by 


itl 


s, 


m 


■£ 


n 


rn^m^+T 


m 


and  precipi- 


(4.8) 


The  distribution  function  of  cloud  droplets  we  denote  as  f  = 
may  express  the  rate  NCp  as  m 


N 


cp 


fk  fm*-k  v(rk‘ 


r 

m*- 


k 


) 


nm/nc ‘ 


We 


(4.9) 


where  v(r-|,  ^)  is  the  collision  kernel  discribed  in  section  4.3. 


The  corresponding  rate  h'cp  is 


m. 


flcp  =  U  <k+m.-k)  fk  Vk  v(rk'Vk  ’ 

(4.10) 


which  may  be  expressed  as 


H  =  (P0/P  )m*N 
cp  °°  *  cp 


(4.11) 


A  realistic  calculation  of  the  cloud  to  precipitation  rate  does  require  some 
information  about  the  cloud  droplet  distribution  function  fm  since  the  size 
m*  lies  in  the  tail  region  of  the  distribution. 


There  is  one  of  the  very  few  exact  solutions  to  the  collisional 
coalescence  problem  which  provides  a  frame  work  for  parameterization  of  the 
CP  process.  The  solution  for  fm  (^  beginning  with  a  single  size  of  cloud 
droplets  r-i  at  time  t=t0  subject  to  a  constant  collision  kernel  v0  is 
(Reference  39) 

fm(t)  =  (1-T)  Tm_1  (4.13) 

where  T  is  determined  by  the  solution  of 


£  -  MO  (1-T)2  (4.13) 

with  wo  =  v0  where  nco  is  the  number  density  of  cloud  droplets  at 

time  t=t  when  collisions  become  to  be  more  important  than  condensation 
growth  in°determining  cloud  droplet  size.  Substituting  the  form  (4.12)  into 
the  CP  rate  expression  (4.9)  we  find 


n2m  v 
C  *  o 


(1-T)2  Tm*~2 


(4.14) 
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we  shall  use  this  result  as  the#basis  for  this  simplified  two  group  model  by 
adopting  the  forms  for  Ncp  and  Hcp  given  by  Eqs.  (4.11)  and  (4.14)  with 
the  following  provisions: 

(1)  The  choice  of  vQ  becomes  an  effective  parameter 
of  the  model-  Loosely  speaking,  it  should  be  se¬ 
lected  as  an  "average"  collision  kernel  over  the 
range  1  <  m  <  m* 

(2)  For  purposes  of  describing  the  conversion  to  pre¬ 
cipitation  rate,  the  cloud  droplet  spectrum  is 
approximated  as  originating  as  a  single  size  at 
the  radius  f-j  =  Rc 


4.4.4  The  Cloud  Collection  (CC)  Process 


The  colli  si onal  interaction  of  precipitation  drops  with  cloud  droplets, 
we  assume,  results  only  in  coalescences  which  enter  the  precipitation  group. 


N 


cc 

We  approximate  this  result  as 
N 


cc 


c  p 


Ncc 

is  then 

formally  defined  as 

ra*  °° 

ncnp 

EE 

m=l  k=m. 

fm  fk  v(rk>  rra) 

,+l 
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v(Rc 

■  V 

(4.15) 


4.4.5.  The  Precipitation  Evaporation  Process 

The  evaporation  of  precipitation  drops  in  unsaturated  cloud-free  air 
has  to  be  considered  to  complete  the  processes  that  balance  the  liquid 
water  existing  as  cloud  droplets  and  as  precipitation  drops.  The  process 
is  represented  by  the  precipitation  evaporation  rate  ftpe-  We  assume 
that  all  precipitation  drops  evaporate  at  a  rate  given  by  that  of  a 
droplet  at  the  average  precipitation  size  Rp. 

Let  us  now  determine  ti.e  rate  hpe  •  The  evaporation  rate  from  a  drop¬ 
let  of  radius  Rp  in  stagnant  air  may  be  expressed  as 

(v)  '  4ljRp°  np  (pvrpvoo)/pco  (4'16) 
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where  D  is  the  diffusion  coefficient,  pVi 

at  the  surface  of  the  particle  and  PVoo  'is 
ment.  v 


is  the  saturation  vapor  density 
the  vapor  density  of  the  environ- 


u.F°r  a  in  a  convective  flow,  the  evaporation  rate  is  enhanced. 

Inri  LC!rre latlons  have  been  developed  to  model  the  increased  evaporation 

and  they  are  generally  of  the  form 


(”pe)  '  (Hpe) 

'  '  conv  '  ' 


(1+C 


stag 


S+> 
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(4.17) 


1 dsnumber  Re  is  based  upon  the  droplet  diameter  and  sedimen- 
tation  veiocity  relative  to  the  air,  and  the  Schmidt  number,  Sc,  is  for 

r -6nV?7fi r  ^1ffusi°n  in.t0  the  ambient.  The  coefficient  C  has  the  value 
C-  0.276  recommended  in  Reference  40.  The  precipitation  evaporation  rate 
is  then  expressed  as 


H_  =  4TrnR„np/S 
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[  H|(H-Hp)  ]  (4.18) 


4.5  Illustration  of  the  CP  Model  for  the  Homogeneous  Cloud 

Vertical  inhomogeneity  is. an  imDortant  aspect  of  any  cloud  and  precipita¬ 
tion  process.  However,  preliminary  to  incorporation  with  the  general,  turbu¬ 
lent,  vertically  inhomogeneous  PBL  model  we  may  examine  some  of  the  character¬ 
istics. of  the  CP  model  for  a  homogeneous  cloud  with  given  turbulence  and  liquid 
Wo l6 r  "inputs* 


We  thus  consider  a  homogeneous  cloud  which  at  time  t  =  0  consists  of  a 
given  amount  of  liquid  water  H^  existing  completely  as  cloud  water 
[Hp(t  -  0)=0]  with  droplets  of  radius  Kc.  We  consider  the  presence  of 
a  uniform  precipitation  flux  divergence  term  which  we  represent  as 


9Xi 


(H  V  ,)  =  H  V  /£ 
p  pi  P  pz  C 


(4.19) 


where  is  a  modeling  parameter  for  this  homogeneous  illustration  only  re 
presenting  an  equivalent  characteristic  vertical  gradient.  We  further  as¬ 
sume  that  the  liquid  water  total  H#,  existing  initially  is  not  replenished 
by  further  decreases  in  saturation  mixing  ratio  as  water  precipitates  from 


To  carry  out  specific  calculations  we  must  specify  the  values  of  the 
collision  efficiency  functions  Els  E2,  E3  and  the  average  collision  kernel 
v0  in  the  cloud-precipitation  conversion  rate,  Eq.  4.14.  We  represent  v 
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as  the  geometric  mean  of  the  extreme  kernels  over  the  collision  range  of 
partners  from  r-|  to  r*: 


v 

o 


ri)v(rl5 


o 


(4.20) 


This  representation  has  the  property  vo  ■+■  o  as  v(r-i,r-i)  ■+  o  as  indeed 
it  should  for  the  model  in  which  the  original  cloud  is  considered  to  con¬ 
sist  of  single  size  droplets.  For  the  collision  frequencies,  we  assume 
Ez  ~  E3  since  both  involve  differential  setting  velocities.  We  utilize 
the  summary  formulae  thats  in  Reference  41  for  these  collision  efficiencies. 
A  fundamental  gap  is  the  absence  of  data  for  collision  efficiency  Ei  of 
drops  colliding  in  the  presence  of  shear.  Although  this  efficiency  should 
depend  upon  the  relative  Reynolds  number  of  the  colliding  partners  (and 
the  shear  rate),  there  seems  to  be  a  reasonable  validation  of  data  over  a 
wide  range  of  shear  with  a  constant  value  E*  =  0.36  (Reference  41). 

We  adopt  this  value  for  these  illustrations. 

The  evaluations  of  the  cloud  and  precipitation  variables  for  a  range 
of  turbulence  dissipation  rates  and  for  a  nominal  liquid  water  level  for 
the  conditions  of  Table  4.2  are  shown  in  Figures  4.1  through  4.4.  The 
general  trends  as  well  as  detailed  structure  exhibited  are  consistent  with 
model  conditions  of  this  illustration.  We  must  emphasize,  however,  that 
these  results  must  be  viewed  against  the  two  fundamental  conditions  of  the 
illustration:  (1)  Cloud  water  He  and  precipitation  water  Hn  originate  from 
a  fixed  total  liquid  water  content  H&  and  (2)  The  cloud  is  homogenous  with 
a  precipitation  flux  divergence  uniformly  distributed  over  the  cloud  given 
by  Eq.  (4.19).  All  cases  are  terminated  when  the  liquid  water  is  reduced 
to  5%  of  the  original  cloud  water. 

The  most  basic  and  general  trend  in  these  illustrations  is  in  the 
decrease  of  the  time  to  reach  a  maximum  precipitation  drop  number  density 
as  well  as  the  time  for  onset  of  significant  precipitation  flux  as  the 
turbulence  levels  rise.  This  result  is  simply  a  manifestation  of  the  in¬ 
creased  conversion  to  precipitation  rate  Ncp  as  the  turbulence  collision 
rate  increases.  The  second  general  trend  is  in  the  magnitude  of  the  pre¬ 
cipitation  flux  ranging  from  0.08  cm/hr  at  the  lowest  turbulence  level  to 
0.35  cm/hr  at  the  highest  turbulence  level. 

It  should  be  noted  that  at  the  highest  levels  of  turbulence 
(e  =  10  m2  sec3)  (Fig  4.5),  a  significant  precipitation  flux  is  established 
on  the  order  of  several  tens  of  minutes.  Such  turbulence  levels  may  be 
characteristic  of  the  dynamics  within  cumulus  clouds  and  it  is  of  interest 
to  examine  the  evolutions  predicted  here  with  those  of  a  "standard"  (albeit 
unverified)  precipitation  model  for  the  cumulus  cases  (Reference  42). 
Kessler  model  results  are  shown  in  Fig  4.6  for  the  same  illustration  con¬ 
ditions  of  Figs  4. 1-4. 5:  The  time  scales  and  general  evolution  seem  com¬ 
parable  to  the  present  model  for  turbulence  levels  e  >  1  m2  sec3.  It 
should  be  noted,  however,  that  the  Kessler  model  inherently  utilizes 
the  Marshal -Palmer  distribution  function  for  rain  drops  as  an  empirical 
input  characteristic  of  rain  from  cumulus  clouds.  The  CP  model  has  no  such 
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emperical  input  restricting  it  to  such  cumulus  parameterizations.  Thus, 
in  the  absence  of  updraft,  the  high  turbulence  levels  of  the  CP  model 
generate  precipitation  of  high  number  density  and  moderate  size 
(Rd  «  130ym).We  believe  this  is  physically  consistent  with  strong  turbu¬ 
lence  in  the  absence  of  updraft.  With  updraft  present,  the  drops  gene¬ 
rated  by  the  CP  process  would  be  maintained  in  contact  with  the  cloud 
droplets  for  a  longer  duration  before  raining  out  and  hence  grow  to  larger 
size.  Thus,  while  the  Kessler  model  by  virture  of  its  parameterization  is 
incapable  of  discribing  the  stratus  case  in  absence  of  updraft,  we  believe 
the  CP  model  when  integrated  with  fluid  dynamic  mean  motion  including  an 
updraft,  would  predict  rain  drop  sizes  consistent  with  the  Kessler  model 
and  the  Marshal -Palmer  parameterization. 

These  illustrations  indicate  that  the  model  and  its  parameters 
exhibit  results  in  terms  of  time  scales  and  magnitudes  of  precipitation 
sizes  and  fluxes  which  are  consistent  with  those  occuring  naturaly  in 
the  atmosphere.  The  turbulence  levels  which  effect  these  results  are 
typical  of  naturally  occuring  turbulence  levels  in  the  atmosphere. 


Table  4.2 


Conditions  For  Illustration  of  the  CP  Model  For  a  Homogeneous  Cloud 


50  ym 
10  ym 


0.001  kg/kg 
1  km 
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Figure  4.1(a) . 


Evolution  of  precipitation  water  mixing  ratio  Hp  for  a  homogeneous 
cloud  with  all  water  existing  as  cloud  water  Hc  at  time  t  =  0. 
Lowest  turbulence  level  of  these  illustrations,  e  =  .001  m2/sec3. 


Evolution  of  average  precipitation  drop  radius  corresponding 
to  Figure  4.1(a). 


Figure  4.1(b). 
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Figure  4.1(c).  Evolution  of  precipitation  flux  corresponding  to  Figure  4.1(a). 


Evolution  of  precipitation  water  mixing  ratio  Hp.  Turbulence 
level  e  =  0.01  m2/sec3. 


Figure  4.2(a). 


Figure  4 -2(b) . 


Evolution  of  average  precipitation  drop  radius  corresponding 
to  Figure  4.2(a) . 


PP  FLUX  (  CM  PER  HOUR  ) 


.  3t 


€  =  0.01  m2/sec3 


.25 


0 


— i - 1 *- 

3  4  5 

T IME(  HOURS  ) 


8 


Figure  4.2(c).  Evolution  of  precipitation  flux  corresponding  to  Figure  4.2(a). 
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Figure  4.3(a). 


Evolution  of  precipitation  mixing  ratio  Hp  for  turbulence 
level  e  =  0.1  m2/sec3. 
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Figure  4.3(b) 


Evolution  of  precipitation  drop  average  radius  corresponding 
to  Fi gure  4.3(a). 
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Figure  4.3(c).  Evolution  of  precipitation  flux  corresponding  to  Figure  4.3(a) 


Figure  4.4(a).  Evolution  of  precipitation  water  mixing  ratio  Hp  for  turbulence 
level  £=1.0  m2/sec3. 


Evolution  of  precipitation  drop  average  radius  corresponding  to 
Figure  4.4(a) . 


Figure  4 ,4( b ) . 
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Figure  4 .4(c) . 


Evolution  of  precipitation  flux  corresponding  to  Figure  4.4(a). 


Figure  4.5(a). 


Evolution  of  precipitation  water  mixing  ratio  Hp  for 
turbulence  level  e  =  10  m2/sec3. 


Figure  4.5(b).  Evolution  of  precipitation  drop  average  radius  corresponding 
to  Figure  4.5(a) . 
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Figure  4. 5(c) . 


Evolution  of  precipitation  flux  corresponding  to  Figure  4.5(a). 
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Figure  4.6(a).  Evolution  of  the  precipitation  water  mixing  ratio  for  the 
homogeneous  cloud  evolution  conditions  of  Figures  1-4  as 
predicted  by  the  Kessler  model  (Reference  42).  There 
is  no  dependence  upon  turbulence  in  the  Kessler  model. 
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Figure  4.6(b).  Evolution  of  the  average  precipitation  drop  radius  as  predicted 
by  the  Kessler  model  corresponding  to  Figure  4.6(a). 
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5.  A.R.A.P.  MODEL  PROBLEMS 


Model  difficulties  may  be  divided  into  four  areas:  those  involving 
the  basic  turbulent  transport  model,  those  involving  additional  physical 
mechanisms  not  incorporated  within  the  model,  those  involving  the  numerical 
complexities  of  faithfully  following  the  modeled  equations,  and  those 
involved  in  obtaining  adequate  field  data  to  rigorously  test  the  accuracy 
of  the  model  and  to  drive  it  under  cases  of  practical  interest.  We 
are  actively  involved  in  all  but  the  last  experimental  area.  We  believe 
we  can  best  contribute  to  this  last  area  by  exercising  the  model  to  determine 
which  parameters  have  the  strongest  influence  on  the  model  behavior  and 
in  this  way  help  to  give  some  guidance  to  observers  as  to  which  parameters 
most  critically  need  to  be  measured.  A  discussion  of  the  first  three 
areas  are  given  in  the  following  sections. 

5.1  Turbulent  Transport  Model 

An  important  problem  area  of  the  turbulent  transport  model  is  its 
inability  to  correctly  predict  the  horizontal  variance  under  conditions 
which  lead  to  a  strong  disparity  between  the  horizontal  and  vertical 
velocity  variance.  As  pointed  out  in  Section  3.1  in  Reference  18,  the 
model  provides  for  a  Monin-Obukhov  similarity  for  the  horizontal  velocity 
which  does  not  exist  in  the  data.  Similarly  the  horizontal  velocity 
variance  predicted  by  the  model  for  free  convection  departs  significantly 
from  the  data  near  the  top  and  the  bottom  of  the  layer  as  shown  in 
Figure  3.3.2  of  Reference  18.  Perhaps  even  more  important  the  model  cannot 
represent  the  transition  to  strongly  stratified  conditions  where  the  flow 
becomes  essentially  two-dimensional.  As  explained  in  Reference  18,  we 
believe  the  principal  difficulty  is  the  anisotropic  nature  of  the 
turbulence  scales.  In  all  of  these  example  problem  areas  the  vertical 
scale  becomes  much  less  than  the  horizontal  scale  in  regions  of  the  flow. 
Lewellen  and  Sandri  (Reference  16)  attempted  a  modification  to  our  basic 
single  scale  model  which  permitted  the  horizontal  macro-length  scale 
to  be  different  from  the  vertical  macro-length  scale.  The  vertical  scale 
was  identified  with  the  single  scale  as  previously  modeled  and  the 
horizontal  scale  was  allowed  to  vary  with  the  mixed-layer  height.  This 
preliminary  attempt  at  a  2-scale  model  did  successfully  eliminate  the 
Monin-Obukhov  scaling  of  the  horizontal  velocity  variance  and  replace  it 
with  a  dependence  on  the  mixed-layer  thickness  which  is  consistent  with  data. 

However,  further  tests  of  the  2-scale  model  as  formulated  by  Lewellen 
and  Sandri  (Reference  16)  have  shown  that  it  lacks  the  generality  desired 
for  adaptation  into  the  basic  model.  The  most  noticeable  failure  was 
that  in  strongly  stratified  flow  it  did  not  force  the  ratio  of  the 
vertical  velocity  variance  to  the  horizontal  velocity  variance  to  approach 
zero  to  permit  a  transition  to  two-dimensional  turbulence.  The  coupling 
between  the  two  time  scales  given  in  the  Lewellen  and  Sandri  (Reference  16) 
formulation  did  not  allow  this  decoupling.  We  are  currently  reformulating 
the  two  scale  jnodel  to  remove  this  problem. 
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If  the  basic  Rotta  term  for  return  to  isotropy  is  reintrepreted  as  a 
tendency  for  axial  symmetry  about  the  3  coordinate  axes,  then  this  part  of  the 
pressure  strain  term  may  be  written  as 


(5.1) 


37  ^  -T-  (6ij  -  S2iS2j)] 


(5.2) 


As  long  as  there  is  the  single  time  scale  for  the  three  return-to-axial 
symmetry  terms,  Eqs.  (5.1)  and  (5.2)  are  identical.  However,  when  one  of  the 
time  scales  is  assumed  to  be  unequal  to  the  other  two,  then  Eq.  (5.2) 
provides  the  basis  for  a  correction  to  the  Rotta  terms,  namely. 


1 

7 


(5.3) 


The  subscript  (1)  represents  the  single  coordinate  direction  about  which  the 
flow  is  assumed  to  be  biased.  The  time  scale  x  is  generally  taken  to  be  the 
eddy  turnover  time  A/q.  The  analogous  turnover  time  in  the  biased  direction 
would  appear  to  be  A1/(uiui) 1/2 .  When  «  t,  Eq.  (5.3)  suggests  that 
energy  will  be  taken  out  of  the  u^ui  components  on  a  time  scale  of 
A,/(u.Ui  )1/2.  This  will  provide  a  tendency  towards  equal ibrating  the  two  time 
scales  even  when  the  length  scales  are  quite  unequal.  We  are  currently 
proceeding  to  rework  the  analysis  of  Lewellen  and  Sandri  (Reference  16) 
following  the  formulation  suggested  by  Eq.  (5.3).  This  approach  appears 
to  have  the  potential  for  duplicating  those  results  and  providing  the 
transition  to  nearly  two-dimensional  turbulence  under  conditions  of 
strong  stratification. 
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5.2  Uncertain  Physical  Processes 

We  have  attempted  to  incorporate  into  the  model  those  processes  which 
appeared  most  pertinent,  but  a  number  of  unmodeled  phenomena  may  play  a  role 
under  some  conditions.  Certainly  the  micro-scale  cloud  physics  addressed  in 
Chapter  4  plays  an  important  role  in  determining  the  rate  of  precipitation  to 
the  surface  under  conditions  conducive  to  rain  or  drizzle.  We  believe  the 
formulation  described  in  Chapter  4  forms  the  basis  for  incorporating  this 
phenomena  into  the  basic  model,  but  at  the  present  time  the  interaction  of 
turbulence  and  cloud  microphysics  remains  an  area  of  considerable  uncertainty. 

Surface  breaking  along  the  air-sea  interface  is  a  process  not  included 
within  the  model  which  is  important  in  determining  aerosol  concentrations  in 
the  surface  layers  and  may  also  be  important  in  determining  the  effective  zo 
of  the  surface  of  the  ocean. 

The  role  of  internal  waves  within  the  atmospheric  marine  boundary  layer 
is  another  uncertain  phenomenon.  They  almost  certainly  exist  a  significant 
fraction  of  the  time.  How  important  their  interaction  with  the  turbulent 
transport  within  the  layer  is,  remains  to  be  determined,  although  the  sample 
calculation  of  Appendix  A  demonstrates  that  there  can  be  a  strong  interaction 
under  some  conditions. 

A  more  speculative  possibility  is  for  a  double  diffusive  instability  to 
exist  along  the  top  of  the  boundary  layer  when  the  air  in  it  is  both  cooler 
and  more  moist  than  the  air  above  it.  The  boost  that  radiation  provides  for 
the  transport  of  heat  provides  a  difference  in  effective  diffusi vities  for 
heat  and  moisture  which  might  drive  a  double  diffusive  instability.  That  is, 
a  finger  of  warm  dry  air  protruding  into  the  boundary  layer  from  above  has  the 
potential  for  cooling  faster  than  the  moisture  mixes  into  it  and,  thus, 
continuing  to  fall  because  it  becomes  heavier  than  the  surrounding  boundary 
layer  air. 

Under  some  unstable  conditions  the  boundary  layer  turbulence  penetrates 
well  into  the  troposphere  in  the  form  of  cumulus  clouds.  These  cumulus  clouds 
play  an  important  role  in  the  coupling  of  a  global  atmospheric  model  to  the 
surface.  This  penetrative  turbulence  is  beyond  the  capability  of  the  current 
model . 

5.3  Numerical  Difficulties 

The  entire  spectrum  of  turbulent  flow  phenomena  might  be  interpreted  as  a 
numerical  problem  since  the  Navier-Stokes  equations  provide  a  rigorous 
formulation.  The  difficulty  is  that,  although  these  equations  may  be  solved 
rather  straightforwardly  with  today's  computing  resources  for  simple  boundary 
conditions,  the  large  range  of  scales  involved  in  the  motions  of  the 
atmospheric  boundary  layer  prohibits  a  complete  detailed  solution.  Since 
these  motions  range  from  the  small  dissipative  eddies,  as  small  as  10-3 
meters,  it  is  essential  that  some  averaging  and  approximating  is  required  for 
any  practical  model. 
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The  present  model  requires  at  least  an  averaging  over  one  spatial 
coordinate  to  form  a  2-D  unsteady  model.  Most  of  our  calculations  have  used 
the  1-D  unsteady  model,  which  goes  even  further  and  averages  over  both 
horizontal  directions.  This  program  which  assumes  horizontal  homogeneity  is 
appropriate  for  many,  if  not  most  conditions  over  the  open  ocean,  provided  the 
role  of  2  and  3-D  motions  are  adequately  parameterized  within  the  turbulence 
model.  But  other  conditions  such  as  that  occurring  along  irregular  coast 
lines  cannot  be  represented  by  either  the  1-D  or  2-D  models. 

Chapter  3  details  our  current  attempt  at  providing  increased  averaging  in 
the  vertical  direction  so  that  we  can  gain  the  potential  for  resolving  some 
completely  3-D  features  with  the  model. 
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6.  FUTURE  PLANS 


The  dynamics  of  the  temperature  inversion  at  the  top  of  the  boundary 
layer  represents  one  of  the  regions  involving  the  most  uncertainties  in 
current  planetary  boundary  layer  models.  Below  this  altitude,  dynamics  are 
dominated  by  turbulence,  while  above  it,  internal  gravity  waves  dominate.  The 
mutual  feedback  between  turbulence  and  waves  in  this  transition  region 
provides  a  complex  interaction  which  controls  such  important  phenomena  as  peak 
values  of  C&.  We  plan  to  continue  the  type  computations  given  in  Appendix  A 
which  allow  the  large  essentially  2-D  features,  either  internal  waves  or 
turbulent  eddies,  to  be  resolved  as  part  of  the  ensemble  mean  motion  while  the 
smaller  scale  turbulence  is  treated  by  the  turbulence  closure  model.  By 
allowing  the  computation  to  resolve  a  significant  part  of  the  interaction, 
dependency  on  uncertainties  in  the  closure  model  are  reduced.  The  results  can 
then  be  time  averaged  to  compare  with  simpler  1-D  models  or  scaling 
relationships.  Principal  effort,  will  be  expended  on  the  2-D  calculations 
but  we  expect  to  complete  sufficient  analysis  of  these  results  to  suggest 
means  for  improving  the  simpler  representation. 

We  plan  to  continue  our  effort  described  in  Section  4  to  incorporate  more 
cloud  physics  into  our  model.  The  current  model  is  based  on  quite  simple 
physics.  Thermodynamic  equilibrium  is  assumed  to  exist  between  liquid  and  gas 
phase  water  at  all  times.  The  liquid  which  exists  is  assumed  to  be  in  the 
form  of  small  droplets  of  specified  size.  The  droplet  size  is  constant  in 
space  and  time.  Only  the  number  density  varies  as  the  liquid  water  content 
varies.  In  actuality,  we  would  expect  the  droplet  size  distribution  to  be 
controlled  by  a  complex  interaction  between  the  turbulent  fluctuations  in 
relative  humidity  and  the  ambient  concentration  of  condensation  nuclei.  Since 
the  droplet  size  distribution  is  such  an  important  variable  in  determining 
visibility  within  a  fog  of  given  liquid  water  content,  we  are  attempting  to 
make  use  of  the  analytical  and  experimental  studies  performed  by  other  NASC 
contractors  on  fog  droplet  dynamics  to  incorporate,  at  least,  some  droplet 
growth  dynamics  within  our  model.  The  important  interaction  between  thermal 
radiation  and  droplet  size  will  also  be  included. 

A  third  task  for  the  future  calls  for  us  to  exercise  the  model  to 
investigate  conditions  which  lead  to  the  initiation  of  fog  or  low-level 
stratus.  Concurrently,  with  extending  the  model's  capability,  we  wish  to 
utilize  it  to  exemplify  phenomena  of  interest  to  the  Navy.  At  least  some  of 
these  calculations  will  be  in  support  of  the  fog  model  evaluation  study  being 
carried  out  by  Cal  span. 

We  also  plan  to  explore  the  feasibility  of  2  major  extensions  of  the 
model.  We  would  not  expect  to  complete  these  extensions  within  the  next  year, 
but  do  expect  to  determine  the  relative  attractiveness  of  proceeding  with 
these  extensions.  The  first  of  these  calls  for  exploring  the  use  of 
second-order  closure  techniques  to  improve  the  parameterization  of  the  effects 
of  cumulus  clouds  in  global  atmospheric  models.  The  intent  would  be  to  view 
the  occurrence  of  cumulus  as  a  penetration  of  boundary  layer  turbulence  well 
into  the  troposphere.  A  first  step  in  this  exploration  would  be  to  become 
more  familiar  with  how  the  effects  of  cumulus  are  parameterized  in  current 
models.  The  final  task  calls  for  exploring  the  extension  of  our  model  to  3 
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dimensions.  This  should  be  an  attractive  alternative  if  the  1-D  integral 
model  described  in  Chapter  3  works  sufficiently  accurately  to  allow  the 
vertical  distributions  of  the  primary  variables  to  be  represented  by  just  a 
few  grid  points.  A  code  which  is  no  larger  than  our  present  2-D,  unsteady 
model  could  then  be  constructed  to  examine  completely  3-D,  unsteady  events 
which  happen  in  real  coastal  environments. 

The  development  and  testing  of  such  a  3-D  Code  based  upon  a  general  •* 

hybrid  treatment  of  the  full  3-D  problem  described  by  second  order  closure 
turbulence  theory  is  a  task  of  major  proportion.  We  feel  that  validation 
of  the  hybrid  method  for  the  general  1-D  homogeneous  case  including  the 
presence  of  an  inversion,  stratus  cloud,  and  radiative  transport  should 
be  demonstrated  before  commitment  is  made  to  develop  such  a  3-D  hybrid/ 
integral  code.  The  development  and  validation  of  the  inversion  layer 
and  moisture  transport  including  the  presence  of  stratus  cloud  for  the 
1-D  homogeneous  case  as  noted  on  page  20  is  our  immediate  near  term  goal 
for  the  hybrid  procedure. 

We  expect  to  pursue  this  research  in  close  cooperation  with  NEPRF 
personnel  in  order  to  be  responsive  to  particular  questions  which  may 
arise  during  the  course  of  this  investigation. 
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APPENDIX  A 


A  Numerical  Study  of  Breaking 
Kel v in-Helmholtz  Billows  Using  a 
Reynolds-Stress  Turbulence  Closure  Model* 

by 

R.  I.  Sykes  and  W.  S.  Lewellen 
Aeronautical  Research  Associates  of  Princeton,  Inc. 
50  Washington  Road,  P.0.  Box  2229 
Princeton,  New  Jersey  08540 


ABSTRACT 

A  two-dimensional  numerical  study  of  breaking  Kel vin-Helmholtz  billows  is 
presented.  The  turbulent  breaking  process  is  modeled  using  second-order 
closure  methods  to  describe  the  small-scale  turbulence,  whilst  the  large  scale 
billow  itself  is  calculated  explicitly  as  a  two-dimensional  flow.  The 
numerical  results  give  detailed  predictions  of  turbulence  levels  and  time 
scales,  and  are  consistent  with  laboratory  and  atmospheric  observations.  Two 
general  predictions  of  the  model  are  that  the  structure  of  turbulent 
temperature  fluctuations  is  very  different  from  that  of  the  velocity 
fluctuations,  the  former  being  much  more  striated,  and  secondly  that  the 
timescale  of  the  growth  and  breaking  process  is  virtually  completely 
determined  by  the  initial  velocity  shear. 

♦Submitted  for  publication  in  J.  Atmos.  Sci.,  September  1981. 
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1.  Introduction 


Shear  instabilities  are  now  widely  recognized  as  a  significant  mechanism 
for  producing  turbulence  and  mixing  in  stably-stratified  fluids.  There  are 
numerous  observations  of  "billow  turbulence"  in  the  atmosphere  and  oceans  (see 
Maxworthy  and  Browand,  1975)  and  these  events  have  been  associated  with 
Kel vi n-Hel mhol tz  instabilities  of  the  wind  profile.  Laboratory  studies  oy 
Thorpe  (1973)  show  similar  vortex  roll-up  features  and  the  generation  of  a 
turbulent  layer  of  fluid  which  mixes  and  spreads  the  initial  shear  layer 
before  decaying  back  to  a  quiescent  flow. 

The  initial  stages  of  the  instability  are  now  well  understood.  Finite 
amplitude  numerical  computations  by  Patnaik,  et  al .  (1976),  Peltier  ,  et 
al.  (1978)  have  confirmed  the  linear  stability  predictions  of  growth  rates  and 
mode  structure  of  the  growing  waves,  and  have  gone  on  to  calculate  the  roll-up 
of  the  vortex  layer.  However,  these  laminar  calculations  at  moderate  Reynolds 
numbers  have  been  unable  to  identify  the  secondary  instability  which  results 
in  the  breakdown  of  the  vortex  into  a  turbulent  layer.  Davis  and  Peltier 
(1980)  made  calculations  up  to  a  Reynolds  number  of  500  based  on  shear  layer 
thickness  but  failed  to  produce  a  secondary  instability.  In  the  latter  paper, 
the  authors  calculated  local  Rayleigh  numbers  in  the  flow,  and  showed  that 
with  a  Reynolds  number  of  500,  a  region  of  very  high  Rayleigh  number  was 
produced  in  the  rolled-up  vortex.  They  speculate  that  this  region  is  in  fact 
convectively  unstable,  but  that  the  most  unstable  modes  will  be  longitudinal 
rolls  since  the  flow  is  strongly  sheared.  This  hypothesis  accounts  for  the 
fact  that  the  instabilities  have  not  been  triggered  in  the  two-dimensional 
numerical  models. 
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The  aim  of  the  present  study  is  to  model  the  breakdown  of  the  rolled-up 
vortex  layer  and  the  subsequent  spread  and  decay  of  the  turbulence  using  a 
second-order  closure  scheme  to  describe  the  three-dimensional,  small  scale 
turbulence  field.  Second-order  closure  models  have  been  used  extensively  to 
study  horizontal ly-homogeneous  atmospheric  boundary  layer  flows,  and  in 
particular  are  capable  of  representing  the  major  dynamical  processes  involved 
in  turbulent  convection  (see  Wyngaard  and  Cot6,  1974;  Lewellen  and  Teske, 
1975;  Yamada  and  Mel  lor,  1975).  [There  are  difficulties  associated  with  the 
turbulent  transport  of  turbulent  kinetic  energy  (see  Zeman  and  Lumley,  1976; 
Andre,  et  al . ,  1976),  but  this  is  not  a  serious  problem  for  the  prediction  of 
the  mean  profiles  of  the  major  second-order  quantities.]  Since  the  model 
calculates  the  growth  of  turbulent  correlations  in  an  unstable  environment 
moderately  accurately  (because  the  growth  is  driven  directly  by  production 
terms  in  the  equations  without  recourse  to  empirical  closure),  we  may 
reasonably  expect  the  results  of  the  calculation  to  have  some  validity.  One 
would  not  have  the  same  expectation  of,  say,  a  mixing-length  closure  model, 
since  the  timescale  of  the  growth  of  the  turbulence  is  the  same  as  the  mean 
flow  timescale,  while  the  mixing-length  model  implicitly  assumes  that  the 
turbulence  is  in  local  equilibrium.  Thus  turbulent  stresses  and  heat  fluxes 
would  be  generated  as  soon  as  the  density  profile  became  unstable  if  the 
mixing-length  were  related  to  the  stability  as  in  the  model  of  Orlanski  and 
Ross  (1973). 

As  we  shall  demonstrate,  the  second-order  closure  model  produces  results 
which  are  certainly  qualitatively  correct.  Quantitative  data  on  turbulence 
quantities  is  not  sufficiently  detailed  to  make  any  meaningful  assessment  of 
model  accuracy,  but  there  are  more  general  quantities  associated  with  the 
breakdown  which  can  be  used  for  comparison.  There  are  estimates  of  the 


timescale  of  the  turbulence  decay  and  measurements  of  the  final  states  from 
the  experiments  by  Thorpe  (1973),  for  example.  The  final  state  of  the  shear 
layer  is  an  important  quantity  in  deciding  the  importance  of  Kel vin-Helmholtz 
instability  as  a  mixing  process  in  stratified  fluids,  since  this  is  a  measure 
of  the  total  amount  of  mixing  produced  by  the  event.  Another  source  of  data 
is  high-resolution  radar  probing  of  the  atmosphere,  where  vortex  roll-up 
events  have  been  reported  by  many  investigators,  as  discussed  in  Chapter  III 
of  the  main  report.  However,  although  it  is  possible  to  estimate  the 
intensity  of  the  turbulent  fluctuations  from  the  intensity  of  the  radar  echo, 
the  measurements  of  the  initial  wind  and  temperature  profiles  immediately 
prior  to  the  instability  are  generally  inadequate  for  a  crucial  test. 

In  spite  of  the  lack  of  detailed  information,  there  is  certainly 
sufficient  data  to  determine  whether  the  turbulence  model  is  reasonably 
accurate,  and  hopefully  some  of  the  numerical  results  will  prompt  more 
detailed  measurements. 


2.  Model  Equations  and  Numerical  Solution 

The  second-order  closure  scheme  employed  in  this  study  has  been  described 
in  detail  by  Lewellen  (1977)  and  we  therefore  present  only  the  final 
equations.  We  consider  incompressible,  Boussinesq  flow,  and  write  the 
Reynolds-averaged  equations  with  turbulent  correlations  denoted  by  an  overbar. 
For  two-dimensional  motion  in  the  (x,z)-plane,  with  z  in  the  upward  vertical 
direction  we  define  a  stream-function  v,  such  that  u  =  av/az,  W  =  -  a^/ax. 
Then  the  vorticity  s  =  V2^,  and  the  equations  of  motion  are 
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In  the  above,  (U,W)  is  the  mean  velocity,  T  is  the  mean  temperature, 
(u,v,w)  is  the  turbulent  velocity  fluctuation,  6  is  the  temperature 
fluctuation,  and  q  =  [u2  +  v2  +  w2]1/2.  a  is  a  length  scale,  and  appears  in 
the  modelled  terms  in  the  Reynolds-averaged  turbulent  correlation  equations. 
The  values  of  the  empirical  constants  are  as  follows:  b  =  1/8,  A  =  0.75, 
Vc  =  0.3,  and  s  =  1.8.  These  values  were  previously  chosen  to  match  model 
results  with  a  set  of  experiments  with  simple  geometries  as  described  by 
Lewellen  (1977). 

For  the  Kel vin-Helmholtz  instability  calculation,  we  begin  with  a  profile 
of  temperature  and  velocity,  viz. 


T  =  T0  +  AT  tanh 


(11) 


U  =  /iU  tanh  4 
6 


(12) 
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Note  that  velocity  and  temperature  both  have  the  same  profile  shape.  The 
initial  conditions  define  a  Richardson  number  Ri  =  (g/T0) (aT6/aU2) ,  which  will 
be  the  basic  dimensionless  number  in  our  experiments.  The  problem  also 
contains  a  Reynolds  number,  but  since  we  have  omitted  explicit  laminar 
diffusive  terms  from  the  equations  of  motion,  we  are  effectively  studying  the 
high  Reynolds  number  limit. 

In  order  to  begin  the  integration,  we  have  to  specify  initial  values  for 
the  turbulence  energy,  q2,  and  the  length  scale  A;  in  the  calculations 
reported  here,  all  other  turbulence  quantities  were  initially  set  to  zero.  A 
perturbation  vorticity  was  also  added  so  that  the  instability  would  amplify 
this  disturbance  and  produce  the  vortex  roll-up.  The  specification  of  these 
initial  conditions  will  be  discussed  in  the  next  section. 

Periodic  boundary  conditions  were  employed  in  the  x-direction,  and  the 
length  of  the  domain  was  chosen  somewhat  larger  than  the  wavelength  of  the 
fastest  growing  linear  mode.  This  choice  is  suggested  by  the  nonlinear 
calculations  of  Patnaik,  et  al .  (1976),  which  show  that  the  longer  waves  grow 
to  larger  amplitude,  and  may  therefore  dominate  the  finite  amplitude 
development  of  the  wave.  For  most  of  our  integrations,  a  domain  length  of  15g 
was  used,  although  some  runs  were  made  with  different  lengths  to  investigate 
the  sensitivity. 

The  boundaries  in  the  z-direction  were  placed  sufficiently  far  from  the 
shear  layer  for  them  to  have  negligible  effect  on  the  flow  development.  A 
zero  normal  gradient  condition  was  specified  on  all  flow  quantities  on  the 
upper  and  lower  boundaries. 

The  equations  of  motion  are  discretized  on  a  finite-difference  grid  which 
is  uniformly  spaced  in  the  horizontal  but  has  finer  vertical  resolution  in  the 
vicinity  of  the  shear  layer.  The  equations  are  solved  using  centered  spatial 
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derivatives,  and  an  ADI  scheme  for  the  temporal  derivative.  The  advection 
operator  conserves  first  and  second  moments  apart  from  time-dependence  errors 

arising  from  the  use  of  explicit  advective  velocity  components  in  the  ADI 

scheme.  The  time-stepping  is  made  more  stable  by  the  introduction  of  a 
coupling  between  the  mean  variables  and  the  turbulent  fluxes  by  the  artifice 
of  an  additional  implicit  eddy  diffusive  flux  which  is  explicitly  subtracted 
when  the  correct  turbulent  flux  is  added  (Lewellen  and  Sheng,  1980).  Provided 
the  time-step  is  sufficiently  small,  the  errors  introduced  by  this  procedure 
will  be  negligible. 

Finally,  the  Poisson  equation  for  the  streamf unction  is  solved  directly 
using  the  decomposition  method  of  Swarztrauber  and  Sweet  (1975). 

3.  Results 

We  consider  an  initial  state,  as  described  in  the  Section  2,  given  by 

T  =  T0  +  aT  tanh  —  (13a) 

U  =  aU  tanh  —  (13b) 

6 

which  defines  a  Richardson  number  Ri  =  (g6AT/T0AU2) .  This  is  in  fact  the 
minimum  gradient  Richardson  number,  defined  by  Rig  =  C(g/T0) (aT/az)]/(aU/3z)2, 
which  occurs  at  z  =  0.  Linear  inviscid  theory  (Miles,  1961)  shows  that  the 
flow  is  stable  to  small  perturbations  provided  Ri  >  0.25. 

We  first  present  detailed  results  from  a  calculation  with  Ri  =  0.1.  In 
this  example,  we  applied  a  sinusoidal  (in  x)  variation  to  the  vorticity  in 
order  to  initiate  the  instability.  The  perturbation  vorticity  amplitude  of 
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O.l/iU/6  is  quite  large;  this  is  because  we  are  mainly  interested  in  the 
finite-amplitude  breaking  of  the  billow  and  the  subsequent  turbulent 
development,  so  that  we  do  not  wish  to  spend  time  computing  the  initial  growth 
from  a  very  small  perturbation.  Checks  were  made  to  ensure  that  the  initial 
perturbation  was  sufficiently  small  that  it  did  not  influence  the 

finite-amplitude  development. 

The  turbulence  energy  q2  was  initially  set  at  10-3aU2  throughout  the 
domain,  with  the  turbulence  length  scale,  A,  set  equal  to  0.46.  The 

turbulence  energy  is  initially  small  enough  to  be  effectively  zero,  but  the 
initial  length  scale  does  have  some  effect  on  the  results.  Variations  in 
these  initial  values  will  be  discussed  in  more  detail  in  Section  6  below. 

The  integration  was  performed  using  a  domain  of  length  156  and  a  height 
of  126  centered  on  the  mid-point  of  the  shear  layer.  The  computational  mesh 
consisted  of  41  x  61  points,  with  41  points  spaced  uniformly  in  the 

horizontal,  and  a  non-uniform  vertical  mesh  giving  a  grid  spacing  of  0 . 1 6  in 
the  central  region  and  expanding  out  to  roughly  0.26. 

Figure  1  shows  a  sequence  of  contour  plots  of  the  dimensionless 

temperature  field,  (T  -  T0)/aT,  at  dimensionless  time  x  =  1.5,  2.8,  4.3,  5.8, 
7.3,  and  11.8,  where  t  =  t(gAT/T0AU).  Corresponding  plots  of  the  small-scale 
turbulent  energy  q2/AU2  are  presented  in  Figure  2.  The  billow  develops 
initially  as  an  effectively  laminar  flow  with  a  growing  vortex  core  in  the 
center,  and  thin  braids  forming  between  the  cores.  The  initial  turbulent 
perturbation  in  q2  is  amplified  slowly  in  the  braids  where  the  shear  is 
highest,  and  the  local  Richardson  number  is  lowest.  At  x  =  2.8  the 
instability  has  vertically  reached  its  maximum  amplitude,  with  temperature 
contours  quite  convoluted  in  the  vortex  core,  and  the  braids  have  become  very 
thin.  At  this  stage  the  temperature  field  is  statically  unstable  in  the  core, 
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and  the  q2  contours  show  growth  of  turbulence  associated  with  the  buoyant 
instability,  although  the  turbulence  is  still  relatively  small  everywhere  with 
the  buoyantly  produced  turbulence  in  the  core  comparable  with  the  turbulence 
levels  in  the  braids. 

At  the  next  time,  t  =  4.3,  most  of  the  temperature  structure  in  the  core 
has  been  mixed  by  the  turbulence  which  has  reached  a  relatively  high  level 
around  the  outer  part  of  the  core;  the  maximum  value  of  q2  is  O.ISaII2.  The 
core  has  already  begun  to  spread  horizontally,  although  the  braids  are  still 
quite  distinct.  At  T  =  5.8,  the  braids  have  been  mixed,  and  the  turbulence 
has  spread  throughout  the  layer  in  the  horizontal.  There  is  evidence  of 
smaller  scale  disturbances  in  the  temperature  field,  and  a  movie  of  the  time 
evolution  shows  them  to  be  travelling  wave-like  disturbances  moving  along  the 
top  and  bottom  of  the  layer,  and  also  weakly  rolling  up  in  a  manner  similar  to 
the  original  instability.  However,  they  are  also  being  mixed  by  the  high 
background  turbulence  so  that  they  do  not  persist  long.  There  is  also  some 
remnant  of  the  original  large  scale  vortex  which  continues  to  slowly  turn  the 
isotherms.  We  can  still  identify  the  regions  of  production  of  the  turbulent 
energy,  since  the  level  of  q2  in  the  center  is  about  twice  the  level  at  the 
edges  of  the  domain. 

The  trend  toward  flattening  the  isotherms,  and  spreading  the  turbulence 
horizontally  continues  until  at  the  latest  time,  x  =  11.8,  there  is  virtually 
no  further  mean  motion,  and  the  turbulence  is  almost  homogeneous  in  the 
horizontal.  The  turbulence  level  has  dropped  to  0.022aU2  at  this  stage,  i.e., 
almost  a  factor  of  10  below  its  maximum  amplitude  which  occurred  around 
x  =  4.3. 
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The  evolution  described  above  is  qualitatively  similar  to  the  turbulent 
breakdown  of  Kel vin-Helmholtz  billows  generated  in  the  laboratory  as  described 
by  Thorpe  (1973).  In  the  tank  experiments,  the  small-scale  turbulence  is 
generated  in  the  vortex  cores,  and  spreads  horizontally  to  amalgamate  the 
billows  at  dimensionless  time  of  roughly  t  =  t0  +  2,  where  t0  is  the  time  at 
which  the  fine-scale  turbulence  first  appears.  Since  t0  ~  3  in  the  numerical 
integration,  the  obsesrvation  that  the  turbulence  merges  horizontally  at  some 
time  between  t  =  4.3  and  t  =  5.8  is  quite  consistent  with  the  experimental 
data.  Furthermore,  the  thickness  of  the  region  of  turbulent  fluid  increases 
slowly  after  merging  in  the  horizontal  and  seems  to  reach  an  equilibrium  level 
near  the  end  of  the  integration.  This  is  again  consistent  with  the 

experimental  observations  of  Thorpe  (1973),  which  show  the  layer  increasing  in 
thickness  up  to  t  =  i0  +  7,  i.e.,  x  =  10  in  our  case.  Also  the 

non-dimensional  height,  RL  =  g(AT.h)/(2AU2)  where  h  is  the  thickness  of  the 
turbulent  layer,  reaches  a  value  of  about  0.4  in  accordance  with  the 

experiments. 

Some  insight  into  the  dynamics  is  obtainable  from  examination  of  the 
energy  budgets.  In  Figure  3,  we  plot  the  roll  energy,  EK,  and  the  small-scale 
turbulence  energy,  EQ,  defined  by 

{[u-u(z)]2  +  w2}dxdz  (14) 


q2dxdz  (15) 
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where 


L 

U  dx  (16) 

o 

We  can  see  the  rapid  initial  build-up  of  roll  energy  as  the  growing 
billow  extracts  kinetic  energy  from  the  mean  shear,  whilst  the  turbulence 
energy  EQ  remains  virtually  constant  at  its  initial  perturbation  level.  At 
roughly  x  =  2.5,  EQ  begins  to  increase  as  the  billow  starts  to  break,  and  EK 
levels  off  at  its  maximum  value.  The  kinetic  energy  is  effectively 
transferred  via  potential  energy  to  the  small-scale  turbulence  which  reaches  a 
peak  at  x  =  5.  The  large  scale  motion  loses  its  energy  very  rapidly  after 
breaking,  while  the  small-scale  turbulence  persists  throughout  the 
integration,  and  exhibits  a  relatively  slow  rate  of  decay.  Thorpe  (1973) 
shows  that  the  turbulence  persists  until  roughly  x  =  x0  +  12  =  15,  at  which 
point  there  appears  to  be  a  collapse  of  the  turbulent  eddies  into  striated 
structures.  The  present  second-order  closure  model  of  the  small-scale  eddies 
provides  a  poor  representation  of  this  final  stage  of  the  decay  of  turbulence 
in  a  stably-stratified  medium,  and  we  have  therefore  not  attempted  to 
integrate  beyond  x  =  13. 

Profiles  of  the  initial  and  final  (x  =  11.8)  local  gradient  Richardson 
number, 


Ri  G 


g  aT/az 
T0  (au/az)2 


(17) 


are  shown  in  Figure  4.  The  final  value  of  Ri q  in  the  layer  is  not  quite 
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constant,  but  lies  between  0.3  and  0.35.  This  is  again  consistent  with 
Thorpe's  estimate  of  the  final  Richardson  number.  This  value  is  significantly 
larger  than  the  critical  value  of  0.25,  so  that  the  mixing-process  apparently 
extracts  sufficient  energy  from  the  initially  unstable  flow  to  mix  the  mean 
profiles  further  than  the  point  at  which  no  more  energy  can  be  extracted.  The 
agreement  between  the  calculation  and  the  experiment  on  this  number  is  some 
confirmation  that  the  model  describes  the  energy  transfers  in  the  initial 
stages  reasonably  accurately.  The  basic  energy  source  is  kinetic  energy 
produced  by  the  billow  acting  on  the  mean  shear,  and  it  seems  likely  that  the 
final  amount  of  mixing  will  be  determined  by  the  amount  of  energy  that  the 
billow  can  extract  before  it  breaks.  If  the  turbulence  model  predicted  the 
wrong  time-scales  for  the  breaking  process,  we  would  be  unlikely  to  obtain  the 
correct  final  state. 

The  evolution  from  an  initial  Richardson  number  of  0.2  is  illustrated  in 
Figures  5  and  6.  A  mesh  of  41  x  61  grid  points  was  used,  as  in  the  previous 
case,  with  a  domain  size  of  156  in  the  horizontal  again,  but  only  86  in  the 
vertical.  Contours  of  temperature  and  small-scale  turbulent  kinetic  energy  at 
x  =  3.0,  6.0,  9.0,  14.6,  and  20.6.  Clearly,  the  billow  is  much  more 
restricted  in  the  vertical  due  to  the  fact  that  there  is  relatively  less 
kinetic  energy  available  in  the  mean  shear  to  drive  the  instability  against 
the  restoring  buoyancy  forces.  The  higher  initial  Richardson  number  also 
makes  the  turbulent  kinetic  energy  production  in  the  braids  less  important,  so 
that  at  x  =  3  when  the  isotherms  are  only  just  being  overturned,  the  maximum 
energy  is  already  in  the  vortex  core.  Furthermore,  the  restricted  size  of  the 
vortex  core  appears  to  result  in  the  turbulence  being  spread  throughout  the 
region  rather  than  being  confined  to  the  outer  half  of  the  vortex  as  in  the 
Ri  =  0.1  case.  The  velocity  in  the  core  is  also  smaller  in  the  Ri  =  0.2  case. 
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so  there  is  less  tendency  for  the  turbulence  to  be  swept  around  the  outer  part 
of  the  vortex. 

As  in  the  Ri  =  0.1  case,  the  turbulent  vortex  core  spreads  horizontally, 
with  small  scale  waves  running  along  the  top  and  bottom  of  the  layer  producing 
local  turbulence  maxima.  Finally,  the  mean  motion  is  totally  suppressed,  and 
the  isotherms  become  horizontal,  and  we  have  a  homogeneous  layer  of 
turbulence. 

The  evolution  of  the  large-  and  small-scale  kinetic  energies  is  shown  in 
Figure  7.  Comparing  this  with  Figure  3,  we  see  that  the  large-scale  roll-up 
extracts  much  less  energy  than  the  Ri  =  0.1  case,  but  the  maximum  in 
small-scale  energy  is  almost  as  large  as  the  large-scale.  Since  the 
small-scale  energy  arises  from  the  potential  energy  created  by  the  roll-up 
process,  this  demonstrates  that  there  is  more  potential  energy  per  unit 
kinetic  energy  in  the  Ri  =  0.2  case.  This  should  not  be  a  surprise,  since  the 
Richardson  number  is  a  measure  of  the  ratio  of  potential  to  kinetic  energy  in 
the  initial  profile,  and  a  linearized  disturbance  will  contain  energy  in  the 
same  proportion. 

We  may  also  note  that  if  the  time  scale  for  the  Ri  =  0.2  case  is  divided 
by  two,  then  a  number  of  features  match  up  with  the  Ri  =  0.1  case.  Defining 
t1  =  x/Ri ,  then  the  large  scale  energy  maximizes  at  x'  ~  50  in  both  cases. 
Furthermore,  the  small-scale  energy  is  down  to  about  half  its  maximum  value  at 
x‘  =  100.  It  seems  that  the  timescale  is  determined  by  the  shear  velocity  and 
the  layer  thickness,  while  the  buoyancy  timescale  determines  the  details  of 
the  roll-up  such  as  core  size  and  energy  partitioning. 

The  profile  of  gradient  Richardson  number  at  the  final  time,  t  =  20.6. 
is  shown  in  Figure  8  for  the  case  with  initial  Ri  =  0.2.  The  final  state  is  a 
region  of  virtually  constant  Richardson  number  at  a  value  of  roughly  0.45, 
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slightly  larger  than  the  case  with  Ri  =  0.1.  The  laboratory  experiments  of 
Thorpe  (1973)  indicate  that  the  final  Richardson  number  is  independent  of  the 
initial  value,  the  final  value  being  between  0.26  and  0.385,  but  this  result 
is  for  initial  Richardson  number  less  than  0.14. 

In  considering  atmospheric  observations,  we  must  recognize  the  fact  that 
the  high-resolution  radars  are  measuring  small-scale  refractive  index  changes. 
The  latter  are  very  closely  related  to  humidity  and  temperature  fluctuations 
rather  than  the  turbulence  energy  itself.  Since  humidity  is  a  scalar,  it 
satisfies  the  same  equation  as  temperature,  and  it  can  be  shown  that  the  mean 
square  humidity  fluctuations  are  identical  to  the  mean  square  temperature 
fluctuation  provided  the  initial  profiles  are  identical.  It  is  therefore 
instructive  to  look  at  the  behavior  of  temperature  fluctuations,  remembering 
that  asymmetries  can  be  introduced  by  having  an  asymmetric  humidity  profile. 
In  fact,  the  length  scale  of  the  energy-containing  turbulent  eddies  also 
affects  the  refractive  index  fluctuations  at  the  short  wavelengths  visible  to 
the  radar,  but  variations  in  the  length  scale  do  not  produce  any  significant 
changes  in  the  pattern. 

Contours  of  the  mean  square  temperature  fluctuation,  q2/aT2,  for  the  case 
with  initial  Ri  =  0.1  are  shown  in  Figure  9  at  t  =  1.5,  2.8,  4.3,  and  8.8. 
The  most  obvious  feature  is  the  highly  localized  nature  of  02.  Initially,  "t? 
is  concentrated  in  the  developing  braids,  where  the  temperature  gradient  is 
very  large.  As  the  billow  begins  to  break,  there  are  local  patches  of  high 
temperature  variance  in  the  mixing  regions  around  the  edges  of  the  core,  but 
these  patches  do  not  spread  throughout  the  vortex  core  as  it  breaks.  Rather, 
as  the  vortex  breaks  down,  the  temperature  variance  vanishes,  and  only  remains 
in  the  braids  at  t  =  4.3.  As  the  braids  are  eroded  by  the  turbulent  mixing, 
the  temperature  variance  decays,  but  remains  in  narrow  bands  which  are  drawn 
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out  horizontally  until  they  are  almost  flat.  This  qualitative  description 
also  applies  to  the  Ri  =  0.2  case. 

The  difference  between  the  temperature  variance  and  velocity  variance  is 
also  evident  by  comparing  Figure  10  with  Figures  3  and  7.  The  integrated 
value  of  the  temperature  variance  tends  to  both  grow  and  decay  more  rapidly 
then  the  integrated  velocity  variance. 

The  contrast  in  behavior  between  the  temperature  variance  and  turbulent 
kinetic  energy  may  be  somewhat  better  understood  by  examining  the  production 
terms  in  the  second-order  turbulence  correlations.  We  have 

_  aUi  2g  _ 

P(q2)  =  -  2uiUk -  +  —  wo  (18) 

axk  'o 

P(wo)  =  -  w2  —  -  uw  —  -  wo  —  -  ue  —  +  -f-  e2  (19) 

9z  ax  az  ax  T0 

P(o2)  =  -  wo  “  -  uo  —  (20) 

az  ax  v  ‘ 

where  P(  )  denotes  production  terms.  The  only  other  terms  in  the  relevant 
equations  are  diffusion  and  dissipation  terms.  Now,  as  the  vortex  core  rolls 
up,  large  gradients  of  mean  quantities  are  generated,  and  consequently  the 
production  terms  are  also  large. 

Once  the  vortex  breaks  the  mean  gradients  are  rapidly  mixed  away,  and 
therefore  there  is  only  diffusion  and  dissipation  of  "o7  which  must  begin  to 
decay  immediately.  However,  there  is  still  production  of  wo  due  to  the  term 
involving  the  temperature  variance.  Thus  heat  flux  does  not  decay  as  quickly 
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as  IP-*  Furthermore,  wb  is  the  main  production  term  in  the  q^  equation  within 
the  core,  so  the  turbulent  kinetic  energy  continues  being  produced  while  e2  is 
decaying.  Thus  q2  grows  in  the  core,  and  eventually  spreads  outward 
horizontally,  whilst  b2  vanishes  in  the  core,  and  only  remains  significant  in 
the  braids,  where  there  is  gradient  production.  When  the  braids  are  eroded, 
the  temperature  variance  is  slowly  diffused  and  dissipated,  while  the  mean 
shear  stretches  the  regions  of  high  into  nearly  horizontal  stratifications. 

The  differences  between  the  q2  and  “b2  patterns  are  thus  quite  plausible. 
We  note  that  this  dynamical  effect  could  not  be  reproduced  by  an  effective 
viscosity  model  which  has  no  equivalent  interaction  between  second-order 
quantities.  We  should  also  note  the  implication  that  high  resolution  radar 
returns  from  atmospheric  billow  events  of  this  type  should  always  show 
highly-structured  patterns  with  narrow  bands  of  high  intensity,  rather  than  a 
homogeneous  layer  of  turbulence.  This  is  quite  consistent  with  the 
observations  reported  in  the  literature.  Figure  11  showing  the  evolution  of 
K-H  wave  breaking  as  received  by  radar  (Browning  and  Watkins,  1970)  shows  the 
same  features  as  exhibited  in  Figure  9. 

4.  Sensitivity  Studies 

There  are  a  number  of  parameters  in  the  numerical  integrations  presented 
in  the  previous  section  which  must  be  chosen  independently  from  the  physical 
constants  of  the  problem.  These  parameters  include  numerical  discretization 
quantities  such  as  the  grid  size  and  initial  conditions  for  the  dynamic 
variables.  We  do  not  include  the  empirical  constants  in  the  second-order 
turbulence  closure  model,  since  these  are  taken  to  be  fixed  from  comparison 
with  other  experimental  data.  It  is  important  to  obtain  some  indication  of 
the  dependence  of  the  results  or  these  external  parameters  if  we  are  to  gain  a 
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useful  understanding  from  the  numerical  integrations.  The  few  main  quantities 
which  need  to  be  chosen  are: 

4.1  numerical  resolution 

4.2  initial  turbulence  level 

4.3  initial  turbulence  length  scale 

4.4  length  of  the  integration  domain 

We  shall  investigate  variations  in  each  of  these  quantities  independently  for 
the  case  with  Ri  =  0.2. 

4.1  Numerical  resolution.  The  integration  presented  in  the  previous 

section  was  carried  out  with  a  mesh  of  41  x  61  points;  this  was  chosen  to 
provide  adequate  resolution  of  the  small  scale  flow  features.  This  fact  is 
demonstrated  by  comparing  the  results  with  those  from  an  integration  using 
32  x  40  grid-points.  The  kinetic  energy  evolution  from  both  integrations  is 
shown  in  Figure  12.  Unfortunately,  the  low  resolution  integration  had  a 
slightly  smaller  initial  vorticity  perturbation  which  is  responsible  for  some 
of  the  differences,  but  notwithstanding  this  initial  discrepancy,  the  overall 
budgets  are  within  about  10%. 

Contour  plots  from  the  two  runs  show  that  the  low  resolution  integration 
has  much  more  numerical  noise  around  the  point  where  the  vortex  is  breaking 
and  gradients  are  highest.  However,  this  does  not  adversely  affect  the 
evolution,  as  can  be  seen  from  the  kinetic  energy  plots,  and  the  two 
integrations  do  remain  very  close  together  throughout  the  period. 
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4.2 


Initial  turbulence  level.  The  remaining  sensitivity  tests  are  all 


carried  out  using  the  low  resolution  model,  so  we  are  comparing  with  the  run 
in  Figure  12(b).  Figure  13  shows  the  kinetic  energy  evolution  using  an 
initial  value  for  q2  of  3  x  10-3  AU2,  which  is  three  times  that  of  the 
previous  run.  The  main  effect  here  is  a  reduction  in  the  large  scale  energy 
by  about  12%,  and  an  increase  of  about  10%  in  the  small  scale  energy  maximum. 
Thus  a  higher  background  turbulence  level  retards  the  development  of  the 
instability  slightly,  and  also  allows  the  vortex  to  break  a  little  more 
quickly,  extracting  more  of  the  potential  energy  from  the  roll-up.  However, 
the  effect  is  not  large,  and  it  can  be  reasonably  assumed  that  the  results  are 
representative  of  small  background  turbulence  levels.  In  fact,  the  initial 
value  for  q2  was  chosen  as  large  as  reasonably  possible  in  order  to  provide 
some  turbulent  mixing  to  prevent  the  braids  from  becoming  too  thin  to  be 
resolved  by  the  numerical  grid.  Since  the  growth  of  q2  in  the  breaking  is  via 
the  buoyancy  term,  it  is  virtually  independent  of  the  initial  value  for  q2, 
due  to  the  rapid  exponential  growth  in  the  initial  stages.  We  therefore  feel 
justified  in  choosing  q2  large  enough  to  control  the  braid  thickness,  as  long 
as  it  is  not  large  enough  to  influence  the  dynamics  of  the  breaking. 

Although  the  higher  value  of  initial  q2,  3  x  10-3  aU2,  appears  to  have  a 
small  effect  on  the  energy  budgets,  there  was  much  less  energy  in  the  short 
numerical  grid  modes,  so  this  run  was  chosen  as  the  base  integration  for 
comparison  of  turbulence  length  scale  and  domain  length  variations. 

4.3  Initial  turbulence  length  scale.  This  is  probably  the  most  difficult 
and  unfortunately  one  of  the  most  sensitive  parameters  to  choose.  The  choice 
of  initial  A  is  not  entirely  arbitrary,  since  it  must  bear  some  relation  to 
the  thickness  of  the  shear  layer,  6,  if  we  are  studying  turbulence  generated 
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by  the  instability  mechanism  itself.  We  chose  the  value  0.46  for  the  runs 
presented  in  the  previous  section;  this  is  a  reasonable  value  for  a 
turbulence  length  scale,  but  cannot  be  justified  too  strongly.  In  order  to 
test  the  sensitivity,  we  carried  out  integrations  with  initial  A  =  0.26  and 
0.86,  and  the  energy  evolutions  from  all  three  runs  are  shown  in  Figure  14. 

Firstly,  the  large  scale  kinetic  energy  is  largely  unaffected  by  the 
choice  of  A.  The  only  difference  is  in  the  oscillation  in  the  decaying  phase 
of  EK,  which  is  accentuated  by  small  values  of  a.  However,  the  small  scale 
energy  does  depend  quite  strongly  on  the  initial  value  of  a.  Changing  A  by  a 
factor  of  2  in  either  direction  from  0.46  changes  the  maximum  turoulence 
energy  by  almost  a  factor  of  2  in  the  same  direction.  The  source  of  the 
problem  is  the  length  scale  equation,  which  does  not  allow  the  length  scale  to 
grow  sufficiently  quickly  during  the  breaking  stage,  so  that  the  turbulence 
does  remember  its  initial  condition  for  a  long  time.  The  length  scale  does 
grow  as  the  billow  breaks,  but  not  so  fast  as  q2  which  is  virtually 
independent  of  its  initial  condition. 

On  this  point,  all  we  are  able  to  say  is  that  0.46  is  a  reasonable 

initial  value  for  A,  and  the  results  of  this  section  illustrate  the  effect  of 
varying  A.  We  should  emphasize  that  the  qualitative  results,  both  in  the 

energy  plots  and  the  details  of  the  contour  fields,  are  not  significantly 

changed  by  this  variation  in  A.  Apparently  the  initial  value  for  A  ought  to 
be  fixed  by  comparison  of  turbulence  levels  with  experimental  observation,  or 
an  improved  dynamical  equation  for  a  is  needed. 

4.4  Integration  domain  length.  We  choose  the  wavelength  of  the  unstable 

mode  we  are  studying  by  fixing  the  length  of  our  periodic  domain.  The  runs  in 
the  previous  section  used  a  length  of  156,  which  was  suggested  by  the  laminar 
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numerical  integrations  of  Patnaik,  et  al.  (1976)  as  giving  the  mode  which  grew 
to  largest  amplitude. 

Figure  15  shows  the  energies  from  three  runs  with  domain  lengths  of  106, 
156,  and  206.  Allowing  for  the  length  of  the  domain,  the  short  domain  clearly 
restricts  the  growth  of  the  billow  significantly;  the  large  scale  energy 
maximum  is  nearly  a  quarter  of  the  156  domain  case  (i.e.,  the  energy  density 
is  about  35%),  and  the  turbulence  energy  is  about  half  (i.e.,  energy  density 
is  75%).  It  is  therefore  unlikely  that  this  mode  would  be  the  dominant  finite 
amplitude  mode,  since  the  longer  modes  would  continue  to  grow.  The  long 
domain,  206,  has  a  significantly  larger  energy  in  the  vortex,  but  allowing  for 
the  length  of  the  domain,  the  energy  density  per  unit  horizontal  length  is 
only  a  little  longer.  The  turbulence  energy  per  unit  length  is  actually 
smaller. 

Thus  the  case  with  a  domain  length  of  156  seems  to  represent  the  most 
efficient  mode  for  producing  kinetic  energy,  although  there  is  not  a  great 
deal  of  difference  with  a  length  of  206. 

5.  Cone! usions 

This  two-dimensional  numerical  study  of  breaking  Kel vin-Helmhol tz  waves 
using  a  second-order  turbulence  closure  model  to  describe  the  small-scale 
turbulence  has  shown  agreement  with  experimental  observations  on  several 
points.  More  importantly,  given  the  lack  of  precise  quantitative  experimental 
data,  the  numerical  integrations  have  suggested  several  general  features  of 
the  dynamics  of  the  fine-scale  turbulence  in  addition  to  the  specific 
quantitative  predictions  of  the  particular  cases  studied.  These  general 
points  await  further  experimental  evidence  to  provide  confirmation  or 
refutation. 
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The  results  demonstrate  the  utility  of  the  closure  model  in  studying  this 
problem;  the  closure  model  contains  most  of  the  important  dynamical  effects 
necessary  for  a  description  of  the  flow.  Furthermore,  it  is  probably  the 
minimum  sophistication  capable  of  describing  the  fine-scale  turbulence,  since 
it  has  shown  that  the  turbulence  develops  on  the  same  time  scale  as  the  mean 
flow,  so  that  equilibrium  assumptions  would  not  be  correct. 

The  numerical  results  indicate  the  correct  qualitative  behavior,  namely 
the  billow  breaks  to  produce  turbulence  in  the  vortex  core,  which  then  spreads 
horizontally  to  form  a  homogeneous  layer.  The  layer  also  spreads  vertically 
to  a  small  extent.  The  time-scale  and  extents  of  the  spreading  are  in  broad 
agreement  with  the  visual  observations  of  Thorpe  (1973).  The  final  value  of 
the  Richardson  number  for  an  initial  Ri  of  0.1  is  consistent  with  Thorpe's 
measured  value  of  0.32  ±  0.06.  However,  Thorpe  notes  no  dependence  on  initial 
Ri  for  Ri  <  0.14,  whilst  our  final  value  for  Ri  =  0.2  is  significantly  higher. 
It  remains  to  be  clarified  whether  the  experiments  would  show  such  a  trend  at 
higher  val ues  of  Ri . 

The  integrations  at  different  initial  Richardson  numbers  suggests  that 
the  advective  timescale  *,/au,  where  %  is  the  wavelength  of  the  billow  and  au 
the  velociity  difference  across  the  shear  layer,  is  the  important  timescale 
for  the  growth,  breaking,  and  subsequent  decay  of  the  turbulence.  The 
buoyancy  timescale  affects  the  type  of  flow  since  the  ratio  of  the  timescales 
is  proportional  to  the  square  root  of  the  Richardson  number.  It  is  quite 
possible  that  the  buoyancy  timescale  is  also  important  in  the  late  stages  of 
the  decay,  when  the  turbulent  eddies  collapse  as  the  mean  density  gradient  is 
re-established.  We  have  not  attempted  to  seriously  study  this  part  of  the 
evolution,  since  the  closure  model  has  deficiencies  in  the  limit  of 
strongly-stratified  turbulence. 
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The  breaking  of  the  billow  itself  is  seen  to  be  generated  by  a  convective 
instability  of  the  overturned  fluid  at  the  top  and  bottom  of  the  vortex  core. 
This  turbulence  is  swept  around  the  vortex,  and  actually  remains  largely  in 
the  outer  part  of  the  vortex  for  a  significant  length  of  time  in  the  small 
Richardson  number  case.  The  turbulence  eventually  spreads  throughout  the 
layer,  and  amalgamates  to  form  a  horizontally  homogeneous  layer. 

Finally,  the  results  indicate  that  the  temperature  variance,  u2.  evolves 
quite  differently  from  the  turbulent-kinetic  energy.  Different  buoyancy 
production  terms  in  the  second-order  turbulence  correlation  equations  imply 
the  rapid  decay  of  "e2  in  the  breaking  vortex  core,  so  that  the  regions  of  high 
remain  highly  localized  in  space.  Effectively,  o2  is  generated  initially 
in  the  braids  of  the  developing  billow,  and  these  long  thin  patches  are  then 
stretched  out  by  the  mean  shear  as  they  diffuse  and  decay  after  the  braids 
have  been  moved  away.  The  evolution  of  the  Kel vin-Helmholtz  wave  breaking  as 
observed  by  radar  in  the  atmosphere  is  quite  consistent  with  the  evolution  of 
the  variance  of  humidity  or  temperature  predicted  by  the  model  calculation. 
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LIST  OF  FIGURES 


Figure  1  -  Isopleths  of  dimensionless  temperature  (T  -  T0)/ATh  for  the  case 

with  Ri  =  0.1.  The  contour  interval  is  0.2,  and  positive 
contours  (i.e.,  lighter  fluid)  are  denoted  by  a  dashed  line, 
(a)  x  =  1.5,  (b)  x  =  2.8,  (c)  x  =  4.3,  (d)  T  «  5.8,  (e)  x  =  7.3, 
and  (f)  x  =  11.8. 

Figure  2  -  Isopleths  of  dimensionless  turbulence  kinetic  energy,  q2/AU2  at 

the  same  times  as  Figure  1,  i.e., (a)  x  =  1.5,  (b)  x  =  2.8, 

(c)  t  =  4.3,  (d)  x  =  5.8,  (e)  t  =  7.3,  and  (f)  x  =  11.8.  Contour 

intervals  are  as  follows:  (a)  0.001,  (b)  0.005,  (c)  0.01, 

(d)  0.01,  (e)  0.007,  and  (f)  0.002. 

Figure  3  -  Evolution  of  total  kinetic  energies,  EK  and  EQ,  for  the  large 

eddy  and  the  small-scale  turbulence  in  the  Ri  =  0.1  case. 

Figure  4  -  Profiles  of  gradient  Richardson  number,  Rig,  for  the  Ri  =  0.1 

case.  Solid  line  is  the  initial  profile  at  x  =  0,  dashed  line  is 
the  profile  at  x  =  11.8. 

Figure  5  -  Isopleths  of  dimensionless  temperature  (T  -  T0)/aT  for  the  case 

with  Ri  =  0.2  at  (a)  T  ■  3,  (b)  T  =  6,  (c)  x  -  9,  (d)  T  «  14.6, 

(e)  x  =  20.6.  Contour  interval  is  0.2. 

Figure  6  -  Isopleths  of  dimensionless  turbulence  kinetic  energy,  q2/AU2,  for 

Ri  =  0.2  at  (a)  x  =  3,  (b)  x  =  6,  (c)  x  =  9,  (d)  x  =  14.6, 

(e)  x  =  20.6.  Contour  intervals  are  (a)  0.0003,  (b)  0.006, 

(c)  0.007,  (d)  0.003,  (e)  0.001. 

F.igure  7  -  Evolution  of  large  scale  and  small  scale  eddy  kinetic  energies 

for  the  case  with  Ri  =  0.2. 
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Figure  8 


Figure  9 


Figure  10 


Figure  11 


Figure  12 


Figure  13 


Figure  14 


Figure  15 


Profile  of  gradient  Richardson  nunber,  Rig,  for  the  Ri  =  0.2 
case,  at  t  =  23. 

Isopleths  of  dimensionless  temperature  variance,  o2/aT2,  for  the 
case  with  Ri  =  0.1.  (a)  x  =  1.5,  (b)  x  =  2.8,  (c)  x  =  4.3, 
(d)  t  =  8.8.  Contour  intervals  are  (a)  0.004,  (b)  0.02, 
(c)  0.02,  (d)  0.002. 

Evolution  of  the  total  volume- integrated  temperature  variance, 
ET,  for  the  case  with  (a)  Ri  =  0.1,  (b)  Ri  =  0.2. 

Schematic  representation  of  the  life  cycle  of  an  individual 
Kel vin-Helmhol tz  billow  based  on  the  data  in  the  earlier  figures. 
Time  progresses  from  right  to  left.  Thick  lines  correspond  to 
the  detectable  clear  air  radar  echo,  which  started  as  a  single 
layer  at  1243  and  finished  as  a  double  layer  at  1258  GMT. 
Schematic  vertical  profiles  of  (Waz)  are  indicated  before  and 
after  the  occurrence  of  Kel vin-Helmhol tz  instability.  (From 
Browning  and  Watkins,  1970). 

Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  (a)  high 
resolution  41  x  61,  (b)  low  resolution  32  x  40. 

Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different 

initial  values  of  q2/AU2,  (a)  10-3,  (5)  3  x  10-3. 

Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different 

initial  values  of  A/6,  (a)  0.2,  (b)  0.4,  (c)  0.8. 

Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different 

domain  lengths,  (a)  106,  (b)  156,  (c)  206. 
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Figure  1  -  Isopleths 
case  with  Ri  =  0.1. 

(i .e. ,  1 iqhter  fluid) 


of  dimensionless  temperature  (T  -  f0)/ATh  for  the 
The  contour  interval  is  0.2,  and  positive  contours 
are  denoted  by  a  dashed  line,  (a)  t  =  1.5. 
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Figure  1(b)  -  x  =  2.8 
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Figure  1(c)  -  x  =  4.3 
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Figure  1(d)  -  x  =  5.8 
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Figure  1(f)  -  t  =  11.8 
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Figure  2  -  Isopleths  of  dimensionless  turbulence  kinetic  energy,  q2/AU2  at 
the  same  times  as  Figure  1,  i.e.,  (a)  x  =  1.5,  contour  intervals  is  0.001. 
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Figure  2(b)  -  t  =  2.8,  contour  interval  is  0.005. 
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Figure  2(c)  -  t  =  4.3,  contour  interval  is  0.01. 
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Figure  2(d)  -  x  =  5.8,  contour  interval  in  0.01. 
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Figure  2(e)  -  x  =  7.3,  contour  interval  is  0.007. 
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Figure  2(f)  -  t  =  11.8,  contour  interval  is  0.002. 
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Figure  3  -  Evolution  of  total  kinetic  energies,  EK  and  EQ,  for  the  large  eddy 
and  the  small-scale  turbulence  in  the  Ri  =  0.1  case. 
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Figure  4  -  Profiles  of  gradient  Richardson  number.  Rig,  for  the  Ri  =  0.1  case. 

Solid  line  is  the  initial  profile  at  t  =  0,  dashed  line  is  the  profile  at  x  =  11.8. 
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Figure  5  -  Isolpleths  of  dimensionless  temperature  (T  -  Tq)/aT  for  the  case  with 
Ri  =  0.2  at  (a)  x  =  3.  Contour  interval  is  0.2. 
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Figure  5(b)  -  t  =  6. 
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Figure  5(c)  -  x  =  9. 
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Figure  5(d)  -  t  =  14.6. 
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Figure  5(e)  -  r  =  20.6 
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Figure  6  -  Isopleths  of  dimensionless  turbulence  kinetic  energy  q2/AU2  Ri  =  0.2 
at  (a)  x  =  3,  contour  interval  is  0.0003. 
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Figure  6(b)  -  t  =  6,  contour  interval  is  0.006. 
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Figure  6(c)  -  t  =  9,  contour  interval  is  0.007 
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Figure  6(e)  -  r  =  20.6,  contour  interval  is  0.001. 


A- 53 


1.0 
0,9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0  A 
0 

0 


Figure  7  -  Evolution  of  large  scale  and  small  scale  eddy  kinetic  energies  for  the 
case  with  Ri  =  0.2. 
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Figure  8  -  Profile  of  gradient  Richardson  number,  RiQ,  for  the  Ri  -  0.2  case,  at  x  20. o 
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Figure  9  -  Isopleths  of  dimensionless  temperature  variance  e^/AT2  for  the  case 
with  Ri  =  1 .  (a)  x  =  1.5,  contour  interval  is  0.004. 
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Figure  9(b)  -  x  =  2.8,  contour  interval  is  0.02. 
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Figure  9(c)  -  t  =  4.3, 
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contour  interval  is  0.02 
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Figure  9(d)  -  x  =  8.8,  contour  interval  is  0.002. 
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Figure  10(a)  -  Evolution  of  the  total  volume-integrated  temperature  variance,  ET 
for  the  case  with  Ri  =  0.1 . 
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Figure  10(b)  -  Ri  =  0.2. 


A 


HEIGHT  km 


DOUBLE 

LAYER 


STRETCHED 

FILAMENTS 


BRAIDED' 

STRUCTURE 


BREAKING 

WAVE 


SINGLE 

LAYER 


1256 


1253 


1246 


I  243  G  M  T 


TIME 


Figure  11  -  Schematic  representation  of  the  life  cycle  of  an  individual 
Kel vin-Helmholtz  billow  based  on  the  data  in  the  earlier  figures.  Time 
progresses  from  right  to  left.  Thick  lines  correspond  to  the  detectable 
clear  air  radar  echo,  which  started  as  a  single  layer  at  1243  and  finished 
as  a  double  layer  at  1258  GMT.  Schematic  vertical  profiles  of  (A6/Az)  are 
indicated  before  and  after  the  occurrence  of  Kel vin-Helmhol tz  instability. 
(From  Browning  and  Watkins,  1970.) 
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Figure  12  -  Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  (a)  high  resolution 
41  x  61. 
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Figure  12  (b)  -  low  resolution  32  x  40. 
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Figure  13  -  Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different  initial 
values  of  q2/AU2,  (a)  10"3. 
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Figure  13(b)  -  3  x  1CT3. 
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Figure  14  -  Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different  initial 
values  of  A/5,  (a)  0.2. 
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Figure  14(b)  -  0.4 


-I— I — L. 


A-69 


Figures  15  -  Evolution  of  eddy  kinetic  energies  for  Ri  =  0.2  with  different 
domain  lenghts,  (a)  105. 
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Figure  15(b)  -  155. 
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Figure  15(c)  -  206. 
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